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1.0 


BACKGROUND 


1.1  Introduction 

In  integrated  Global  Positioning  System/Inertial  Navigation 
System  (GPS/INS)  navigation  systems,  position  and  velocity  errors 
are  typically  estimated  with  respect  to  some  reference  point  within 
the  host  vehicle.  For  example,  this  may  be  the  point  afforded  by 
the  origin  of  the  accelerometer  triad  or  the  center  of  mass  of  the 
host  vehicle.  The  conventional  approach  is  to  relate  measurement 
residuals  to  error  states  under  the  assumption  that  the  antenna 
phase  center  coincides  with  the  reference  point  while  compensating 
for  the  effects  of  any  lever  arm. 

In  reality,  the  antenna  phase  center  is  usually  at  a  point 
different  from  the  reference  point.  Measurement  residuals 
(observations)  formed  by  differencing  predicted  and  measured  pseudo 
range  and  delta  range  must  incorporate  corrections  for  the  antenna- 
phase-center  lever  arm.  Knowing  the  placement  of  the  GPS  antenna 
with  respect  to  the  host  vehicle  and  the  host-vehicle  attitude  with 
respect  to  some  convenient  reference  frame  provides  the  required 
information.  Eut ,  as  the  conventional  observation  matrix  is  based 
on  coincidence  of  the  antenna  phase  center  and  the  reference  point 
(notwithstanding  that  in  reality  it  may  be  offset),  the  pseudo¬ 
range  and  delta-range  measurement  residuals  are  modeled  to  have 
coupling  only  to  position,  velocity,  user  clock  bias  and  drift 
error  states.  In  this  conventional  approach,  potential  coupling  to 
attitude  errors  or  attitude  rate  errors  is  ignored. 

In  the  conventional  approach,  the  needed  attitude  for  antenna- 
lever-arm  correction  is  generally  provided  by  estimation  of 
attitude  errors  and  propagation  of  the  corrected  attitude  "total 
state"  (typically  a  direction  cosine  matrix  or  quaternion)  to  the 
time(s)  of  measurement.  Attitude-error  estimation  for  such  an 
implementation  depends  entirely  on  correlation  between  the  attitude 
error  and  the  velocity  error.  Specifically,  contribution  of  the 
attitude  error  to  the  velocity  error  is  due  to  the  product  of  the 
attitude-error  vector  by  the  skew- symmetric  matrix  of  specific 
forces  within  the  dynamics  F  matrix. 

This  approach  has  been  successfully  applied  to  many 
applications  such  as  terrestrial  navigation.  There  are,  however, 
applications  where  the  correlation  between  attitude  error  and 
velocity  error  vanishes.  For  example,  when  an  exoatmospheric  host 
vehicle  is  free-falling,  the  specific  forces  are  zero  for  all 
practical  purposes.  Then,  although  GPS  may  succeed  in  estimating 
position  and  velocity  errors  fairly  well,  attitude  estimation  is 
not  possible.  Moreover,  the  correction  for  the  antenna  lever  arm 
becomes  flawed  as  attitude  errors  grow  due  to  instrumentation 
errors  such  as  gyro  bias  and  scale  factor  errors. 
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Typically,  GPS-based  attitude  error  estimation  methods  rely  on 
receiving  carriers  from  multiple  GPS  satellites  and  processing  the 
single  difference  of  their  phases  (Van  Graas  &  Braasch,  1991-2; 
Keirleber  &  Maki,  1991;  Satz  et  al .  ,  1991).  A  minimum  of  3 

noncollinear  antennas  are  required  in  most  methods.  The 
implementation  of  these  methods,  however,  requires  considerable 
receiver  hardware  design  features  and/or  special  software 
processing  beyond  the  typical  GPS/INS  designs. 

An  alternative  approach  to  attitude  estimation  is  proposed  in 
this  report.  We  seek  to  exploit  the  attitude  information 
inherently  present  in  the  pseudo-range  and  the  delta-range 
measurements  obtained  with  a  single  offset  GPS  antenna  mounted  on 
a  platform  undergoing  attitude  changes.  The  attitude  information 
is  recovered  with  standard  recursive  estimation  techniques  such  as 
the  Kalman  Filter  algorithms  of  the  GPS/INS  systems.  The  attitude 
information  is  primarily  recovered  through  processing  delta-range 
measurements.  The  pseudo-range  measurement  is  less  useful  because 
of  the  generally  high  code  loop  noise.  The  implementation  of  this 
approach  does  not  require  any  special  receiver  hardware  and  only  a 
few  minor  changes  to  the  standard  GPS/INS  Kalman  Filter  software. 
Specif ically ,  the  required  changes  to  the  Kalman  Filter  software 
are  simply  a  few  additional  entries  in  the  observation  matrix  H. 

The  proposed  approach  is  ideally  suited  for  applications  where 
the  antenna  is  mounted  on  a  spinning  platform  such  as  a  spinning 
vehicle  in  exoatmospheric ,  free-fall  conditions.  Since  this 
approach  imposes  a  minor  implementation  cost,  it  readily  allows  use 
of  existing  systems.  For  "good"  performance,  an  integrated  GPS/INS 
system  is  required;  the  INS  provides  for  propagation  of  the 
attitude  estimates  in  between  GPS  updates.  In  some  specific 
applications,  however,  GPS  alone  may  suffice. 

An  overview  of  the  contents  of  this  report  is  as  follows. 
Sections  2,  3,  and  4  contain  independent  mathematical  developments 
of  this  approach  to  attitude  estimation.  Section  2  provides  a 
general  mathematical  development  of  the  observation  matrix  H,  which 
is  applicable  to  any  changing  host  vehicle  attitude.  The  H  matrix 
is  developed  without  reference  to  any  of  the  other  models  that  are 
part  of  a  GPS/INS  Kalman  Filter.  Section  2  also  provides  a  quick 
comparison  of  this  method  to  the  more  familiar  "interferometric" 
GPS-based  attitude  determination  methods.  Section  3  provides  a 
mathematical  development  of  the  observation  matrix  H,  which  is  also 
generally  applicable  to  any  changing  host-vehicle  attitude.  This 
development,  however,  is  based  on  a  state-space  formulation  and  as 
a  consequence  it  employs  and  thereby  relies  on  prior  developments  - 
specifically,  the  state  transition  matrix  and  the  plant-noise 
model.  Section  4  provides  a  simpler  analysis  of  attitude 
estimation  but  limited  to  a  spinning  body  scenario.  As  a  result  it 
facilitates  gaining  certain  insights  to  and  identifying 
sensitivities  in  spinning  body  applications .  Section  5  provides 
simulation  results  for  a  spinning-body  application  by  using. 
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the  more  general  formulation  of  Section  3.  Section  6  provides  a 
summary  and  conclusions. 

1.2  Notation 

Vectors  in  specific  coordinate  frames  will  generally  be  of 
lower  case  and  always  with  a  single  underline.  Coordinate  free 
forms  will  be  denoted  with  an  over-bar.  Matrices  will  generally  be 
of  upper  case  and  always  with  a  double  underline.  Scalars  may  be 
uppei  or  lower  case  but  without  any  underline.  Normal  superscript 
designation  of  the  reference  frame  will  be  suppressed  when  the 
frame  is  the  preferred  frame.  The  subscript  "n"  will  be  appended 
to  terms  for  the  nominal  trajectory. 
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2.0 


GENERAL  DERIVATION  OF  THE  OBSERVATION  MATRIX  H 


2.1  Introduction 

This  section  provides  a  general  mathematical  development  of 
the  observation  matrix  H,  which  is  applicable  to  any  changing  host- 
vehicle  attitude.  The  pseudo-range  (PR)  and  delta-range  (DR) 
measurement  residuals  obtained  from  the  corresponding  measurements 
and  their  prediction  at  the  antenna  phase  center  are  related  to 
estimation  quantities  referred  to  the  reference  point.  This  is 
done  to  capture  potential  attitude  and  attitude  rate  information 
lost  in  the  more  typical  process  of  computationally  bringing  the 
antenna  phase  center  to  the  reference  point  (but  using  lever-arm 
corrections  in  formation  of  the  residuals) .  If  coupling  to 
attitude  exists,  then  the  attitude  of  a  host  vehicle  undergoing 
attitude  changes  as  well  its  position  and  velocity  can  be  estimated 
even  if  the  specific  forces  are  zero. 

The  following  analysis  first  derives  H  matrix  row  vectors 
corresponding  to  pseudo-range  residuals  for  a  single  "offset"  GPS 
antenna  and  then  proceeds  to  derive  H  matrix  row  vectors  for  delta 
range  residuals.  GPS  pseudo-range  measurement  performance  is  based 
on  the  code  tracking  loop,  while  that  for  the  delta  range  is  based 
on  the  corresponding  time  increments  of  the  carrier  loop.  It  is 
recognized  that  for  limitations  in  the  antenna-le^er-arm  length, 
the  code  loop  noise  will  generally  overwhelm  the  relation  of 
pseudo-range  residuals  to  attitude  error.  The  delta  range 
residuals,  however,  being  based  on  the  very-low-noise  carrier 
tracking  may  well  afford  estimation  of  attitude  errors  of  the  free 
falling  host  vehicle  (as  well  as  provide  estimation  of  attitude 
rate  error  so  that  gyro  errors  may  be  observed)  .  The  derived  H 
matrix  includes  rows  for  both  pseudo-range  and  delta-range 
residuals  for  completeness  and  improvement  that  may  accrue  from 
pseudo-range  measurements  over  a  period  of  time  (reduction  of  the 
effect  of  code  loop  noise  through  integration) . 

The  more  general  H  matrix  derivations  are  simply  extended  to 
multiple  antennas  but  this  is  of  little  value  to  the  problem  of 
interest  as  it  would  entail  considerable  hardware  complication  with 
attendant  cost.  A  brief  investigation  of  higher  order  terms  in  the 
temporal  DR  residuals  process  is  included. 

A  rotating  antenna  by  itself  introduces  (by  virtue  of  its 
rotation)  carrier  phase  shifts  beyond  those  associated  with  the  DR 
measurements  including  a  fundamental  phase  change  due  to  the 
antenna  rotation  (+_  360  degrees  of  phase  shift  per  360  degrees  of 
rotation  in  the  antenna  plane)  and  phase  anomalies  peculiar  to  the 
antenna  phase  pattern.  These  phase  shifts  must  be  identified  and 
compensation  made.  It  is  assumed  that  the  INS  can  measure  the 
antenna  rotation  with  sufficient  accuracy  so  that  the  unwanted 
phase  shifts  as  well  as  other  antenna  phase  imperfections  can  be 
compensated.  The  problems  and  issues  associated  with  the 
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corresponding  compensation  of  delta-range  residuals  and  antenna 
calibrations  were  not  investigated  here. 

2.2  Pseudo-Range  Measurement 

The  geometry  of  the  problem  is  shown  in  figure  1.  The  reader 
should  be  aware  of  the  direction  taken  for  the  unit  vector,  which 
is  shown  from  the  jth  GPS  satellite  to  the  phase  center  of  the 
offset  GPS  antenna  on  the  host  vehicle.  This  may  be  of  opposite 
sense  to  some  definitions  used  by  others.  This  unit  vector  is  also 
well  approximated  by  a  unit  vector  from  the  jth  satellite  toward 
the  user  as  the  length  of  the  antenna  lever  arm  : s  very  small  in 
comparison  to  the  range  to  the  satellite.  Likewise,  in  this 
analysis,  the  sign  of  error  states  is  such  that  when  added  to  the 
predicted  ^propagated)  "nominal"  total  state,  a  new  estimation  of 
the  total  state  is  provided  (which  becomes  the  new  nominal  total 
state  in  an  extended  Kalman  filter) .  The  pseuao  range  between  the 
antenna  phase  center  and  the  jth  satellite  is  given  by 

z3=(^u  +  -  r3j)  •  e:  +  Bu  ,  (1) 

where  "ru,  ry,  and  "rS]  represent  coordinate  free  vector  form.  As  bias 
in  the  satellite  clock  is  not  modeled,  it  is  not  included  in  (1) . 
In  (2)  we  have  coordinatized  with  respect  to  some  convenient 
reference  frame  Xa .  Frame  notations  are  such  that 

Xe  -  Earth  frame,  generally  Earth  Centered  Earth  Fixed  ( ECEF ) 

X;  -  Inertial  frame,  generally  Earth  Centered  Inertial  (ECI) 

Xn  -  Local  level  frame 

Xt  -  Body  (host  vehicle)  frame. 

As  a  random  process,  the  pseudo  range  is  given  by 

z  =  eaT(ra  +  ra  -  ra  )  +  B  +  v_, 

3  — ]  — u  A  — s]  u  ' 

or 

(2) 

=  e}  (ru  +  rA  -  r£j)  +  Bu  +  vFRj 

when  suppressing  the  superscript  "a"  for  notational  clarity.  We 
will  frequently  suppress  the  superscript  "a"  for  our  selected 
convenient  reference  frame  but  maintain  superscripts  for  other 
specially  designated  frames. 
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ry  is  the  random  process  vector  providing  the  location  of  the  GPS 
antenna  but  as  represented  in  Xa.  ry  will  be  rapidly  changing  due 
to  the  changing  attitude  of  the  host  vehicle.  However,  when 
represented  in  the  host  vehicle  frame  (as  rj  A  )  »  it  will  be  constant 
and  related  by 


r  (t)  =  Ca (t)  rb  . 

n  — b  A 


(3) 


It  is  to  be  noted  that  the  direction  cosine  matrix  Ca  is 

modeled  as  a  random  process  but  not  rjy,  which  is  taken  to  be  well 
known  and  constant. 

For  purposes  of  the  derivation  of  the  H  matrix  for  an  extended 
(or  linearized)  Kalman  filter,  we  assume  a  continuous  spatial  as 
well  as  temporal  distribution  of  z3 ,  ry,  and  ry.  In  the  approach  of 
this  section  it  is  convenient  to  disregard  pseudo-range  and  delta- 
range  measurement  noise.  The  spatial  differential  which  we  will 
take  with  respect  to  the  nominal  or  predicted  value  (i.e.  evaluated 
along  the  nominal  trajectory),  corresponds  to  the  pseudo  range 

measurement  residual.  That  is,  8{-)  =  (•)  -  (•)„  where  we  note 


5  (£t£) 


=  £T£-Ea£n  =  (£n  +5£)  T  (£n  +  5£) 

-  ( 5j3 )  T£r  +£nT$£  =  +  • 


In  particular, 
5z 


>z,  = 


5[eT(r  +  r  -r  ,)]  +  5b„ 

— -]  — u  — A  —S3  u 

(r 


r  -r  .  )T6e.  +  eT  5(r  +r  -r  .)  +  5b 

— An  — sjn  — 3  — jn  — u  — A  — sj  u 


(4: 


As  e.j  is  of  constant  (unit)  magnitude,  5^  for  small  magnitudes  must 
be  orthogonal  to  ey  and,  therefore,  to  ry  +  r*  -  to  which  _e. 

is  parallel.  Hence, 

8za  =  eT  8r  +  eT  8r  +  8bu  (5) 

5z •  =  eT  8r  +  eT  (8ca)  rbn  +  8b 

J  — ]  n  — u  —3  n  srsfc,  — An  11 

as  rbn  =  rb  is  a  constant  (column)  vector. 

— An  — A 
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In  the 
that  of  Zb , 


case  of  8ca  where  the  spatial  variability 

b 

8ca  =£b  "£b  =£b  (£tn~i=)  . 


is  due  to 
(6) 


For  small  displacements  of  Ey,  from  Zbn  (the  nominal  body 
frame),  is  skew- symmetric .  We  define 

54/bn  =  84^  ^  =  C*n -I  .  ( 7 ) 

For  notational  convenience,  8  is  used  in  this  definition  but  this 
does  not  allow  for  separation  of  8  from  84^  as  this  would  presume 
that  *F  has  meaning  for  large  displacements  of  2^  from  Ebn .  8¥  is 
not  defined.  We  then  have 

8ca  =  Ca  84/!in  .  (8) 

— -bn^^= 


Therefore, 


-££; 


igipnrbn 


6b,. 


(9) 


But  as 


Qa  glphn^-bn  =  Qa  g\j/bn  q  bn  £  a  ^-bn  _  g\p  n 


(10) 


and  since 


8¥  r^-R  8£, 


R  = 

*An 


h:n 


"rAzn 

0 


Ayn 

-r 


Axn 


-r 


Ayn 


0 


we  have 

8z.=eT  8r  -eT  R  8$  +  8Bu  • 

:  — ]n  — u  — Dnssftr, — *-  u 


(ID 


(12) 


R  is  the  skew- symmetric  matrix  form  of  vector  r  and  80  is  the 

=An  — An  - *~ 

vector  form  of  skew-symmetric  84|  .  Note  that  8  cannot  be  separated 
from  _80_  in  the  context  used. 
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The  H 
to  the  GPS 


matrix  for  pseudo-range  measurements 
antenna  phase  center  is  given  by‘ 


only  but  referred 


8z, 


e!T 
— In 


JL<t>  Sbu 

-e;  R  1 


(13) 


el 

— kn 


-e*  R 

kn^p^ 


If,  instead  of  bringing  the  H  matrix  to  the  antenna  phase 
center,  had  we  derived  an  H  matrix  based  on  the  antenna  phase 
center  coinciding  with  the  reference  point,  that  is  R  =  0  (but  still 

using  measurement  residuals  based  on  antenna-lever-arm 
corrections) ,  the  resulting  H  matrix  would  show  no  coupling  between 
the  residuals  and  .  By  eliminating  this  noncoupled  term,  we 
have 


8z, 

8z2 


8_r 


8bu 

l 

l 


(14) 


8z 


k 


eT 

— kn 


l 


which  is  the  usual  H  matrix  encountered  when  pseudo-range 
measurements  are  presented  to  a  Kalman  filter  implemented  in  ECEF . 


The  more  general  H  matrix  of  (13)  appears  to  have  promise  for 
observability  for  attitude  but  would  require  a  lc  ge  offset 
distance  of  the  antenna  phase  center  from  the  reference  point  of 
estimated  quantities  to  overcome  the  noise  in  the  GPS  receiver  code 
loops.  Common  sense  dictates  that,  with  a  single  antenna, 
regardless  of  the  number  of  GPS  satellites  for  which  pseudo-range 
measurement  are  made  at  the  same  time,  full  determination  would  not 
be  possible  without  prior  information  from  other  measurement  times. 
The  H  matrix  will  not  be  of  full  rank  for  a  single  antenna  scheme. 
In  particular,  three  separate  times  of  measurement  with  the  offset 
antenna  occupying  three  different  noncollinear  positions  are 
required  for  full  attitude  observability,  i.e.,  ensure  full  rank  of 
(13)  .  If  we  constrain  the  attitude  variability  of  the  free-falling 
host  vehicle  to  spin  about  a  known  inertially  stable  axis,  then 


Remember  that  suppression  of  superscripts  for  vectors 
applies  only  to  Z , . 
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only  two  measurement  times  with  noncoincident  antenna  positions 
should  be  required  for  full  observability.  With  multiple  antennas 
(not  considered  here) ,  unique  or  overdetermination  at  a  single 
measurement  time  is  certainly  possible. 

2.3  Delta-Range  Measurements 

With  respect  to  the  pseudo-range  measurement  residuals,  we 
have  a  first  order  approximation  to  a  spatial  difference  process 
(which  involves  difference  with  respect  to  a  nominal  trajectory). 
Proceeding  onward  to  include  delta-range  measurements,  we  now  seek 
the  spatial  difference  process  of  a  temporal  difference  process. 
We  will  maintain  a  first  order  approximation  to  the  spatial 
difference  part  in  order  to  use  a  linear  Kalman  filter.  The 
temporal  difference  will  be  initially  given  a  first  order 
approximation.  Later,  a  second  order  approximation  will  be 
examined  to  see  what  additional  terms  (preferably  involving 
existing  error  states)  have  coupling  to  the  delta  range 
measurements . 

Taking  the  time  differential  of  the  pseudo-range  measurement 
residual  we  have 

d(5zj  =  d(eT  8r  )  -  d(eT  R  5b)  +  d  (  5b  )  ,  (15) 

3  —]n  — u  —  jn=An — u 


and  on  expanding. 


d(5zj)  =  ( )  8r_  +  ejn  d  ( 5 r_u ) 

-  (deT  )  R  5b  -eT  (dR  )  5b  -  er  R  d  (5b)  +  d  (5b,)  . 

—m  =An — *■  —on  =An  ~ ^  — J”=An  — X  J 


(16) 


By  linearity  we  can  exchange  d(*  )  and  5(  -  )  for  all  except  5b  , 
which,  as  previously  noted,  is  not  separable  with  any  proper 

meaning.  For  5£  we  must  return  to  the  defining  relations.  From 
(7)  and  (8),  on  considering  the  temporal  behavior  we  have 


5ca  =  ca  54fn  ,  54fr  =  £bn  -  ^ 

84*  A  84^  =  Ca  (Cbn-I)Cbn  =  CaCbn  -  I 
D(84/)  =  (DCa)  Cbn  +  C a  ( DCbn )  =  ( DC a )  Cb  +  Ca  (DCbn)  Si 

d  { 5T)  =  Qa  ==na  ”aQa  ~h=-  cT  =  Si 

*ab  ^^>n,a  *ab  ^a.bn  — ab  — a,  bn  — ab  * 
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Q  is  the  skew-symmetric  matrix  form  of  angular  velocity  of  the 

host-vehicle  frame  with  respect  to  the  selected  convenient 
reference  frame.  The  corresponding  vector  form  is 


D(5d>)  =  to  -to  =  5to 

- —  — «b  —a ,  bn  — 

(18) 

d(80)  =  StoLdt  , 


which  is  the  mathematical  statement  of  the  fact  that  the  rate  of 
attitude  error  is  equal  to  the  angular  velocity  error. 

It  is  to  be  noted  here  that  5  is  separated  from  G)lt ,  as  angular 
velocity  is  a  vector  quantity  without  further  qualification.  The 
relations  of  (18)  are  intuitively  obvious  but  we  have  been  careful 
to  preserve  formality.  Since  our  focus  will  be  on  resolution  of 
SCO,,  into  gyro  error  states'  it  is  noted  that  if  we  take  I,  to  be 
the  actual  frame  and  not  "what  we  think  it  is"  irrespective  of  our 
selection  for  Za,  then 


8ci)  =  to  -  to  =  (to  +io) 

— ab  — ab  —a,  bn  — -a  i  — if. 


(to  +  to  ) 

-—a  i  — ;  ,  bn 


8co  =to  -  to  =  5(o  . 

““ab  —lb  —  i ,  bn  — ib 


From  (16)  and  (19)  we  now  have 

8(dz  )  =  (deT  )  8r  +  eT  (8v  )  dt 

J  — 3  n  — u  — }  n  — u 

-  (deT  )  R  80  -  eT  (dR  )  80 
-eT  R  (8(0  Jdt  +  8f  dt , 

— ^  rises  at*  ib  u 


(19) 


(20) 


where  we  have  used  d_r 


9r 

-T_-idt  =  v  dt  and  dB.,  = 

OF  — u 


9b 

— ,  d  dt  -  f  dt  . 
c)t 


2  Resolution  into  gyro  error  states  does  not  imply  that  these 
error  states  will  be  observable. 

J  v^  is  the  velocity  of  the  user  relative  to  the  center  of 
Earth  as  determined  by  a  Z,  fixed  observer.  v.,  is  denoted  here  as 
terrestrial  velocity  when  Z,  =  Ze,and  v,.  is  denoted  as  inertial 
velocity  when  Za  =  Z, .  For  Za  =  Zn,  a  local  level  frame,  v,  is  much 
less  useful.  We  proceed,  however,  along  this  path  and  make 
adjustments  later  for  expression  in  terms  of  terrestrial  velocity 
regardless  of  the  selection  for  Za. 
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Using  the  Theorem  of  Mean  Value  for  integration  and 
takmg  5r.  and  5^  to  be  constant  over  the  interval  (t:-T,t.)  ,  the 
DR  measurement  residual  is  given  by 


(dz  )  =  5  f  dz,  =  (  f  deT  \6r  +  eT  f  8v  dt 

J t,-T  ]  \J t,  -  T  ~ >7  — > u  — —J 

-[  P‘  deT  )r  50  -  eT (  f'1  dR  \  50 

\  J t,-T  — 37— An— —3^  Jt  -T  —An/  ~ ^ 

-eT  R  I  ‘  5co,  dt  +  f  ‘  5f„dt  , 

3nssss AnJt  -q1  ^  _q- 


(21) 


where  the  other  terms  not  included  within  the  integration  are 
evaluated  at  some  undetermined  time  within  this  time  interval.  On 
performing  this  integration  as  far  as  we  can  yields 


5  (Az3) 


=  AeT  5r 

— J  n  — -u 


-(AeT  R  +  e7  AR  \50 

\  —]nm «ar  — Jn  —An  /  — 

-eT  R  f  ‘  5(0  dt  +  f  '  8fdt  . 

:n— AnJt  -T  l*>  Jr  - t 


(22) 


We  arrive  at  a  first  order  approximation  to  the  temporal 
difference  part  if  we  consider  the  navigation  reference  point,  the 
local  clock,  and  attitude  changes  to  be  temporally  "well  behaved" 
so  that  8v  ,  5(0  .  and  5f  ,  are  regarded  to  constant  over  (t.-T,  t,)  . 

— u  —  1  b  u  1 


Then 


5(Az  )  =  AeT  5r  +  eT  t5v 

J  — j  n  “"u  —  j  r. 

-  (AeT  R  +  eT  AR  )  50 

—  3™ An  —313  —An  - 

-eT  TR  5co  +  t5  f, 

—  Jn  —An  — it1  J 


(23) 
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or,  more  specifically, 

5 ( Az^  ( t A )  )  =  AeJ,  (t; )  8_ra  (t.)  +  e_!n  (t  J  T8v 

-(AeT  (t,)R  (tb)  +  eT  (t.)  AR  ( t , )  )  5<t>  ( t . )  (24) 

-eMtd)  TRta(t  e)  6co;fc  (t .)  +  T8fu(tJ  , 

where  t;-  T  <_  ta,  tb,  tc,  td,  te  <  t,  . 

Incorporating  both  pseudo-range  and  delta-range  residuals,  the 
resulting  H  matrix  is 


8ru 

8v 

8£ 

8co 

— ib 

8b.; 

Sf, 

5z1 

eT 

—In 

0T 

-e;  R 

1  n  —A J) 

_oT 

l 

0 

8z2 

'T 

— 2n 

0/ 

~e;R 

*-nsAn 

0/ 

l 

0 

8zk 

— kn 

_oT 

-eT  R 

— R«=An 

0T 

l 

0 

8  (Azj) 

AeJ 

ej  T 

—In 

-  (AeJ  R  +  e;  AR  ) 

~ej  TR 

— =An 

0 

T 

8  (Az2) 

AeJ 

_eTt  T 

-(Ae.TR  +  eTAR  ) 

— ^.'nsAn  — *.n  *Ar. 

-eT  TR 

— in 

0 

T 

8  (Azk) 

AeT 

— kn 

eT  T 

— kn 

-(AeTR  +eTAR  ) 

— 1 ^®=An  — 1 =Ar, 

-eT  TR 

— r.n  wtftj, 

0 

T 

If,  instead  of  bringing  the  H  matrix  to  the  antenna  phase  center, 
had  we  derived  an  H  matrix  based  on  the  antenna  phase  center 

coinciding  with  the  reference  point,  that  is  R^=0  (but  still  by 

using  measurement  residuals  based  on  antenna-lever-arm 
corrections)  ,  the  resulting  H  matrix  would  show  no  coupling  between 

the  residuals  8<j)  and  5o)ifc  .  By  eliminating  these  noncoupled  terms, 
we  have 
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5ru 

5v 

5bu 

5fu 

8z, 

'  eT 

—in 

1T 

l 

0 

5z; 

ej,n 

oT 

l 

0 

5zk 

eT 

—ten 

_0‘ 

1 

0 

5  ( Azj ) 

AeT 

—In 

-eT  T 

— In 

0 

T 

5  (Az2) 

AeT 

-eT  T 

— 2n 

0 

T 

5  ( Azk) 

AeT 

— kn 

-el  T 

— kn 

0 

T 

which  is  the  usual  H  matrix  encountered  when  both  pseudo-range  and 
delta-range  measurements  are  presented  to  a  Kalman  filter 
implemented  in  ECEF. 

As  indicated  earlier,  when  Xa  =  Xn  we  generally  want  the 
terrestrial  velocity  of  the  user  rather  than  the  velocity  as 
perceived  by  an  observer  fixed  to  Xn .  In  coordinate  free  form  and 
considering  only  temporal  variation,  we  define  the  terrestrial 
velocity  as 

Deru  =  v  •  (27) 


By  the  theorem  of  Coriolis, 


=  D  r 

+  0)  a 

x  r 

a  u 

ea 

i 

=  D  r  -  co 

e  u  ea 

x  rx 

(28) 


or,  in  matrix  form 


Dr 

— u 

Vu 


=  CMDr*)  - 

=  v  -  IT 


Q  r 

u 


(29) 


Since  we  take  Xa  to  be  the  actual  frame  (X,  does  not  have  a 
nominal  counterpart  in  this  analysis) ,  we  have 

5 v^  =  5v  -  D  5_r  _  (30) 
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(31) 


Using  terrestrial  velocity  v,  (23)  becomes 
8(AzJ  =  ( AgU  -  ej  TQ  )  5^  +  t5v 

-  (AeT  R  +eTAR  )5b  -  eT  TR  8(0  +  T8f,  . 

— —>n  =An  —  — ;n  =An  — -r  J 


To  a  first  order  AeT -eT  TQ  =  (CaAee  );  which  is  the  temporal  change 

— jn  — jn  =ea  =e  — jn  ^ 

m  the  satellite  unit  vector  as  perceived  by  an  Earth  fixed 
observer  and  expressed  in  Sa . 


then  have 

as  defining 

the  observation 

matrix 

H, 

8_r  8v 

8£ 

8co  v 

8bu  Sf,, 

8z, 

c* 

e;n  £T 

-e*nR 

0_T 

1 

0 

02, 

eTn  07 

-eT  R 

— 

r 

1 

0 

8zk 

0T 

-eT  R 

— 1 

oT 

1 

0 

5  ( Azj ) 

(CaAe1en)T  e^T 

-  ( AeJ  R  +eTAR  ) 

-e‘rTR_ 

1 

T 

5  (Az2  ) 

(C^Ae:n)T  eJnT 

-(Ae‘R  +  ej  Ar  ) 

"n*An  — An 

-eT  TR 

l  n  —An 

0 

T 

5  (Azk) 

(C^Ae^)T  e-'nT 

-  ( AeJ  R  +  eT  AR  ) 

— 1 <»An  — i ^  *Ar. 

-e^TR^ 

0 

T 

If  we  model  8co  t  to  consist  only  of  gyro  bias  5co  and  scale 
factor  4  So)t  / 


8(0  =  £  8ob  + 


(33) 


4  Thus  far  we  have  structured  the  error  state  elements  to  be 
partitioned  into  3d  vectors  and  we  will  continue  to  do  so  here  for 
convenience.  One  may  want  to  model  only  the  x  component  of  scale 
factor  or  only  certain  g-sensitive  terms  or  the  like. 
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where  we  note  the  underlining  t  and  jxo  to  he  single  vector 

quantities  and  not  5  (j_)  and  d*-r:::> 


Wb 


(0° 

— lbz 


(34) 


Using 

part 


terrestrial  velocity,  we  then  have  for  the  delta 

5  (Ac  )  =  (Ae_:  -e_:  TQ  )  8 r_  +  ey  T8v 
-<Ae_~R  +  e/tAP.  j  Sp 

-  e  ;t  TR^  C’  5w;  -e“  TR  C’Wr8o);  +  T5f  . 


range 

(35) 


The  H  matrix  of  expanded  column  dimensionality  would  then  be 
constructed  with  these  additional  terms. 

2.4  Higher  Order  Terms 

Up  to  this  point,  we  have  used  only  a  first  order 

approximation  to  the  temporal  difference  process  A(*  ).  A  first 

order  approximation  is,  of  course,  used  for  the  spatial  difference 

process  5(*  )  ,  to  obtain  a  linear  form  for  the  Kalman  filter.  While 

we  retain  the  first  order  approximation  to  the  spatial  difference 
process,  we  here  briefly  investigate  whether  other  coupling  terms 
may  arise  from  a  higher  order  approximation  to  the  temporal 
difference  part. 
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Repeating  (22) ^ , 


8  (Az  )  =  AeT  8r 

J  —  j  n  — u 


+  eT  f  8v  dt 

— 3nJt,-7  ~ u 


+  eT  AR  \86 

— =*=An/ - 


-eT  R  f  '  8co  dt  + 

— jn"AnJc  -t  — lb 


(22) 


In  taking  8v„,  8o>.b,  and  8fu  to  be  constant  over  (t,-T,  t,)  to 
arrive  at  (23)  is  tantamount  to  approximating  each  of  these 
functions  of  time  by  the  first  term  of  a  Taylor  series  expansion 
about  tj.  We  now  look  at  additional  terms. 

Considering  only  temporal  variation,  a  Taylor  series  expansion 
of  8v„  about  t4  and  by  interchanging  8  and  time  differentiation 
yields 


8v  (t)  =  8v  (td  +8 

— u  — u  1 


dv  , 

"3T1 


(t  -t;)  +  h.o.t. 


(36) 


Relating  to  coordinate  free  form, 

0vu  <— >  D;ru  =  Dfru  -  toiax  (cojaxFu)  -  2toiaxDaru 

Da^u  =  De*u  ~  “ae^u  =  (37) 

a  A  Dfru  ,  V  =  Deru  , 

a  is  the  inertial  acceleration  and  v  ,  as  previously  noted, 
is  the  terrestrial  velocity. 


•'  It  will  be  more  convenient  to  initially  use  v,  in  this 
development  rather  than  terrestrial  velocity  v  and  then  make  the 
adjustments.  We  seek  to  arrive  at  an  acceleration  error  state  that 
may  be  related  to  existing  accelerometer  error  state  elements  of 
the  state  vector. 


We  then  have 


8(Dv  )  =8/a  -  Cl  Cl  r  -2Q  (v  -  Q  r  )\ 

— u  \ —  *ia“ia — u  «=ia  —  =ea — u  / 

=5 a.  -  2Q  8v  -  (Cl  -  2Q  Q  )  8_r 


(38) 


By  using  (30),  (36),  and  (38)  we  have 

8  vu  ( t )  =  5  v  ( t , )  -  Q  8_r  u  ( t  j ) 

+  ( t  -t  ■)  [  8a  ( t : )  -2CI  8v  ( t , )  -  (CF  -2CI  Cl  )  8r  ( t . )  ]  , 

1  —  =1  a  *  =1  a  —i  a=ea  u 


(39) 


and 


8vu ( t ) 


-  fiQ  +  (t-td  (Q2  -2CI  Cl  )  ]8r  (t;) 

|_^sa  =*ia  =ia=ea  J  u 

+  [l-2  (t-ti)Q  J8v(t1) 

+  (t-tj  5a  ftj)  . 


(40) 


We  then  have  the  second  order  approximation  for  the  integral 


ri8- 


8v  dt  =  - 


Q  T-  (Q2  -2Q  Q  )  _L 

ea  fcia  “iaa*!ea  2 


5r  (t, ) 


+  ( IT  +  £2  T2)  Sv(tJ  --^_8a  (tj  . 
**  ""ia  “  2  ~ 


(41) 


In  coordinate  free  form. 


a  =  D2ru  «  f  +  G 

or,  in  matrix  form,  the  corresponding  errors  are 

8  <a  =  5_f  +  8  G  , 

0 

where  f_b  is  the  specific  force  vector  which  is  measured  by  the 
accelerometers  and  G  is  the  gravitation  vector.  As  8g  is  not 
modeled  here,  we  delete  this  term.  If  we  take  Sco^  and  Sfu  to  be 
constant  over  (t^T,  tj  ),  from  (22)  and  (41)  we  have 
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5 (Az) i ) 


jAe^  , 


Q  T-  (ft  -2 ft  ft  )  Z_ 

— »..m  — i  a  — l  a— icte  2 


+eT  ( iT+ft  T- )  8v 


-{Ae_Jr|R  +  eJAR 
-eT  R  t8co 

—  K=An  — it) 

-eT Zl8f 

-:r.“2-  - 

+T8f . 


>  80 


(42) 


Though  not  pursued  further  here,  Sf.b  is  expressible  in  terms 
of  accelerometer  error  states  such  as  bias,  scale  factor, 
misalignments,  etc. 

A  similar  approach  could  be  taken  with  regard  to  8co,b,  not 
making  the  assumption  that  it  is  constant  over  ( t x  — T ,  t  j )  .  However, 
an  angular  acceleration  error  state  would  result  that  is  not 
modeled.  Moreover,  if  it  were  modeled,  an  issue  arises  on  the  time 
propagation  of  the  corresponding  total  state  because  we  do  not 
generally  incorporate  corresponding  sensors  that  would  aid  its 
propagation  in  a  high  dynamic  environment. 


2.5  Comparison  To  Standard  Attitude  Measurement  Techniques 

Although  a  comparison  to  an  exhaustive  list  of  other  GPS-based 
attitude  determination  methods  is  beyond  the  scope  of  this  report, 
a  quick  comparison  of  the  method  described  in  this  report  to 
"interferometric"  attitude  determination  methods  (Van  Graas  & 
Braasch,  1991-1992)  may  be  useful  at  this  point.  We  can  construct 
a  configuration  which  affords  a  convenient  general  comparison 
between  the  technique  presented  here  and  the  more  familiar  attitude 
measurement  techniques.  For  purposes  of  comparison  we  take  any  two 
noncoincident  GPS  antennas  of  the  "standard"  interferometric  methcH 
and  the  set  of  their  phase  differences  for  each  of  four  GPS 
satellites  to  comprise  a  spatial  GPS  Interferometer.  A  set  of 
delta-range  measurements  taken  over  the  same  time  interval  ior  four 
GPS  satellites  using  a  single  antenna  in  effect  provides  a  temporal 
version  of  a  GPS  Interferometer  as  defined  above.  While  more 
restrictive  than  need  be  due  to  use  of  recursive  estimation 
techniques  that  preserve  information  by  propagating  it  forward,  the 
typical  implementation  of  a  conventional  five  channel  GPS  receiver 
having  pseudo-range  and  delta-range  measurements  would,  in  fact, 
conform  to  such  a  configuration. 

A  brief  comparison  between  the  technique  described  in  this 
report  (denoted  as  "Delta-Range-Based  method")  and  standard  GPS 
interferometric  methods  is  provided  in  table  1. 
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Table  1.  Comparison  of  standard  interferometric  GPS-based  attitude 
measurement  techniques  to  delta-range-based  method. 


Standard  (Multiple  Antennas) 

DR-Based  (Single  Antenna) 

Requires  multiple  antennas  and 
a  special  GPS  receiver  design. 
No  requirement  on  host -vehicle 
motion  relative  to  GPS 
satellite  constellation. 

Requires  only  a  single  antenna 
and  a  standard  GPS  receiver 
having  pseudo-range  and  delta- 
range  measurement  capability. 
Requires  changing  attitude  of 
host  vehicle  relative  to  GPS 
satellite  constellation. 
Integrated  GPS/INS  generally 
needed  for  good  performance. 

Any  two  non.coi ncident  antennas 
receiving  carriers  of  four 

GPS  satellites  and  the  single 
difference  of  their  phases 
constitute  a  spatial  GPS 

Int  erf erometer . 

Each  set  of  DR  measurements 
for  four  GPS  satellites  using 
a  single  antenna  on  ho°t 
vehicle  undergoing  attitude 
changes  mav  be  viewed  as 
constructing  temporally  a  GPS 
Interferometer  having  as 
antennas  the  start  and  stop 
an senna  locations.  The 
accumulated  phases  of  the  DR 
measurement  set  are  equivalent 
to  the  set  of  differences  of 
received  phases  at  the  start 
and  stop  antenna  locations. 

A  minimum  cwo  nonpar  1  lei 

GPS  Interfere  -.ers  (a  minimum 
of  "-.hree  non^.  linear 
antennas)  is  required  for 
attitude  determination  at  any 
''  -.  ime . 

A  minimum  of  two  nonparallel 
equivalent  GPS  Interferometers 
(two  DR  measurement  sets) 
required  for  attitude 
determination  with  INS 
propagation  of  information 
over  the  DR  time  period. 

Attitude  is  determined  for  end 
of  period  and  is  propagated 
forward  by  INS. 

The  satellite  phase  errors  are 
removed  by  the  single  spatial 
difference  in  phase  between 
the  two  antennas  of  the  GPS 
Inter feromet  er . 

Removal  of  satellite  phase 
errors  is  inherent  to  each  DR 
measurement . 
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Standard  (Multiple  Antennas) 


DR-Based  (Single  Antenna) 


Ambiguity  resolution  for 
integer  cycles  must  be 
accomplished  for  the  single 
phase  differences.  Use  of 
redundancy,  search  procedures, 
triple  differences,  etc.,  are 
required . 

No  ambiguities  in  each  DR 
measurement  provided  no  loss 
of  carrier  lock  or  cycle  slip 
in  GPS  receiver  channel 
associated  with  GPS  satellite. 

May  require  calibration  and 
compensation  for  antenna 
induced  error  due  to  phase 
pattern  differences  between 
the  two  antennas  constituting 
the  GPS  Interferometer. 

Fundamental  phase  change  due 
to  attitude  changes  in  single 
antenna  must  be  compensated 
for.  Antenna  must  also  be 
calibrated  for  phase  pattern 
anomalies . 

Double  difference  (difference 
between  two  independent  single 
differences)  required  to 
remove  receiver  timing  errors, 
electrical  path  bias  errors, 
etc . 

Use  of  single  antenna  and 
standard  DR  measurement 
obviate  the  need  for  double 
differences . 

Direct  solution  of  attitude 
performed . 

Recursive  estimation  (using 
Kalman  Filter)  of  attitude 
errors  relative  to  nominal 
trajectory  established  by  INS 
and  its  corrections  (resets) 
is  performed. 

2  J 


3 . 0 


STATE -SPACE  FORMULATION 


3.1  Introduction 

In  developing  a  measurement  model,  the  method  presented  in 
section  2  has  taken  a  direct  approach.  The  strength  of  such  an 
approach  is  that  it  has  remained  independent  of  any  other  analysis 
and,  thereby,  has  been  exploratory  in  seeking  coupling  between 
pseudo-range  and  delta-range  measurement  residuals  and  the  error 
states.  The  results  for  the  pseudo-range  residuals  need  no  further 
elaboration  and  are  useful  as  derived.  With  regard  to  delta-range 
residuals,  however,  a  weakness  is  that  it  has  not  been  established 
where  certain  terms  in  (22)  need  to  be  evaluated  to  be  equivalent 
to  (21)  .  For  example,  when  the  time  interval  T  over  which  the 
temporal  difference  is  taken  is  small  such  that  the  change  in  R  is 

small,  the  problem  is  minimal  and  any  reasonable  choice  of 
evaluation  time  will  provide  useful  results.  For  greater 
efficiency  in  coupling  to  the  attitude  errors,  however,  we  wish  to 
maximize  Ar  .  The  choice  of  where  to  best  evaluate  R  in  the 
interval  (t^T,^)  is  uncertain. 

In  this  section,  we  utilize  a  presumed  existing  dynamic  model 
to  develop  the  temporal  differences  required  for  the  delta-range 
relations  through  backward  transition  matrices.  The  strength  of 
this  approach  is  that  there  is  no  uncertainty  as  to  where  to 
evaluate  the  terms  in  the  resulting  relations.  The  weakness, 
however,  is  that  it  depends  on  a  previously  developed  dynamics 
model  and  is  subject  to  any  inadequacies  in  this  model.  In 
particular,  uncertainty  in  the  backward  time  propagation 
contributes  to  uncertainty  in  the  measurement  model. 

3.2  Reference  Frames 

The  previous  approach  required  careful  identification  and 
manipulation  of  representation  in  various  reference  frames.  The 
state-space  formulation  requires  only  minimal  consideration.  The 
only  assumption  made  with  regard  to  reference  frames  for  the  error 
state  vector  corresponding  to  the  previously  developed  dynamics 
matrix  is  that  position  errors  and  attitude  errors  are  represented 
in  the  same  selected  preferred  reference  frame.  In  this 
development,  position  vectors  and  the  skew-symmetric  matrix 
corresponding  to  the  antenna-lever-arm  position  vector  will  also  be 
represented  in  this  preferred  reference  frame.  Otherwise,  the 
error  state  vector  can  mix  and  match  reference  frames,  the  only 
requirement  being  consistency  with  the  previously  developed 
dynamics  model. 


3.3 


State  Equations 


The  system  model  will  be  taken  to  be  nonlinear  but  with 
additive  noise  (Maybeck,  1979;  Anderson  &  Moore,  1979).  In 
particular,  the  total  state  dynamic  model  is  taken  to  be  continuous 
and  given  byb 

x(t)  =  f_[x(t),t]  +G(t)w(t)i  (43) 

where  w(t)  is  a  zero-mean  white  Gaussian  noise  process.  For  the 
purpose  of  either  a  linearized  or  extended  Kalman  filter,  we  assume 
a  nominal  trajectory'  of  similar  form 

xn(t)  =  f_[xn(t)  ,  t]  .  (44) 

The  spatial  difference  process 

8x ( t )  =  x ( t )  -  xn ( t )  (45) 

has  the  first  order  approximation 

8x(t)  =  F[t;xn(t)  ]  8x(t) +  G (t)w(t)  ,  (46a) 


where 


F[t;X|(UUai(|;t)  (46b) 

In  the  case  of  an  extended  Kalman  filter  coupling  a  GPS 
receiver  with  an  IMU,  the  nominal  trajectory  is  updated  (reset) 
following  measurement  update  to  the  new  total  state  estimate  so 
that  the  error  state,  now  reset  to  zero,  need  not  be  propagated. 
As  indicated  in  the  previous  development,  propagation  of  the  total 
state  is  largely  accomplished  through  use  of  IMU  data. 

The  discrete  time  measurements  are  taken  to  be  a  nonlinear 
function  of  the  state  with  an  additive  white  noise  sequence. 

z  (tj  =h[x(ti),tj  +v(t1).  (47) 

Associated  with  the  nominal  trajectory  is  the  sequence  of  nominal 
measurements 

2.  n  ( t ; )  =  h  [xn  ( 1 1 )  ,  tj  .  (48) 


b  It  is  noted  that  (43)  implies  definition  of  "attitude  total 
states."  Attitude  total  states  are  not  usually  defined.  The 
formalism  in  (43),  however,  conveniently  leads  to  the  definition  of 
attitude  error  states. 
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The  spatial  difference  process 

5z(tt)  =  2.  (t.)  -z  (ts)  =h[x(t1),t]  -  h  [x  (t,),t;]+v(t)  (49) 


has  the  first  order  approximation 

8  Z_  ( 1 1 )  =  H[ti;xn(ti)]8x(ti)  +v(ti: 


where 


r  ,  v  ,  dh  [x  ( t; )  ,  t ;  ]  . 

H[t1;xn(t1)]  (ti, 


;50) 


If  we  express 

h'  [SxCtJ  ,x(t  ,)  t  J  =  h[x(ti)ti]  ~  h[xn(tj)  ,  t,] , 
where  we  treat  Sx  as  an  independent  variable,  then 


8:<*0 


Further  discussion  and  proof  is  given  in  Appendix  A. 
For  this  problem  we  can  express 


h,  (x,  t4) 

h, '  (8x,x,  t,) 

h(x,  t,)  = 

h;  (X.t,) 

and  h'  (8x,x,  t1 )  = 

h; '  ( 8x,  x,  t . ) 

hm'  (8x,X,  Tj  ) 

so  that  the  jth  row  of  Htt^-x^  (tt)  ]  is  given  by 

,  Tr .  .  .  8h  (x,  tj  9h. '  (8x, x, t, )  , 

*  — Tx— ~~d (Tx~ 


(51) 


(52) 


:53) 


:54; 


3 . 4  Approach 

The  approach  is  to  develop  a  spatial  difference  of  the 
temporal  difference  process  over  the  time  interval  (t_-T,  tj  for 
the  jth  delta  range  (i.e.,  the  delta  range  for  the  jth  satellite). 
The  delta-range  residuals  are  determined  in  terms  of  the  geometry 
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in  figure  1.  We  then  use  a  previously  developed  linearized 
dynamics  model  to  express  error  state  terms  having  time  arguments 
of  tt-T  by  the  product  of  the  backward  transition  matrix 
(appropriately  premultiplied  by  other  matrices  to  isolate  the 
required  terms)  and  the  error  state  vector  having  the  time 
argument  t,.  We  then  have  a  "row"  of  the  measurement  model 
utilizing  the  argument  5x.  We  then  apply  the  results  of  Appendix 
A  and  (54)  to  obtain  the  corresponding  row  of  the  H  matrix  in  the 
perturbation  model  of  (50) . 

3.5  Delta  Range 

From  figure  1  we  have  the  vector  from  the  jth  satellite  to  the 
antenna  phase  center  at  any  time  t  and  represented  in  some 
preferred  frame  Xa  as 

r  (t)  =  r  (t)  +  r.(t)  -  r  (t)  .  (55) 

The  rf  range  Rj(t)  from  the  antenna  phase  center  to  the  jth 
satellite  (or  that  from  the  jth  satellite  to  the  phase  center  as  we 
take  range  to  be  a  positive  value)  must  account  for  the  range 
equivalent  user  clock  offset 

R3  (t)  =  |  r.  (t)  ||  +  Ba(t)  .  (56) 

We  also  have  nominal  rf  range  corresponding  to  the  nominal 
traj  ectory 

Rjn(t)  =  «r.n(t)||  +  Bun(t)  .  (57) 

For  the  rf  range  and  nominal  rf  range  at  time  t-T,  we  would  have 
the  same  forms  as  (56)  and  (57)  but  with  arguments  t-T. 

The  measured  delta  range  as  a  discrete  time  process  (at  t.) 
and  corrupted  by  measurement  noise  is 

zDRj(t,)  =  R,  (t.)  -  Rj  ( t  j -T)  +vrR.(t;).  (58) 


Although  (58)  has  the  appearance  of  the  difference  of  two 
ranges,  it  must  be  remembered  that  the  measurement  is  not  obtained 
as  the  difference  of  two  pseudo-range  measurements  but  rather  as  a 
receiver  measure  of  the  accumulated  carrier  phase  (integrated 
Doppler  shift)  from  t,-T  to  t,.  In  this  regard,  v:F,(t.)  is  the 
discrete  measurement  noise  process  involving  the  effective  carrier 
tracking  loop  noise  over  the  same  interval . 
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The  nominal  delta  range  corresponding  to  the  nominal  trajectory  is 


z:,K;r  ( t ; )  =  Rjn(t:)  -  R,n  ( t ,  -T )  . 


(59) 


The  delta-range  measurement  residual  discrete  process  is  then 

8ZDRj(^"i)  _  ZDRj^i^  ~  ZDRjn^i^  (60) 

=  Rj  ( t , )  -  R:n  ( t , )  -  [Rj  ( t  i  -T)  -R3n(t,-T)]  +vDRj(t;). 

As  for  any  time  t 

R:n(t)  =  flr^U)  +rAn(t)  -r,.  (t)  ||  +  Bun(t)  (61) 

and  as  spatial  difference  processes  using  8  (•)  =  (•)-(•)  n,  are 


8r_  (t)  =  r  (t)  -r  (t) 

5rA(t)  -  r(t)-7>)  (62) 

8Bu(t)  =  Bu(t)-Bun(t)  , 

we  have 

Rjn(t)  =  II  ru  ( t )  -  5ru(t)  +  jrA(t)  -8rA(t)  -  (t)  ||  +  Bu(t)  -5BJt)(63) 

*B.r,  (t)  -  SrMt)  "  8r>A(t)  ||  +  Bu(t)  -8Bu(t). 


Hence  the  delta  range  measurement  residual  discrete  process  for  the 
jth  satellite  is 


5zDR3  (t,)  =||r:  (t,)  ||  -|r .  (t,)  -8ru(t;)  -8rA(t,: 


-{I  r  .  (t  i-T)  II  -  lx,  (t i-T)  _5£u(t;-T)  '^^(tj-T)  ||  +8Bu  ( t .  -T)  | 

+  vrP,3  <  t , )  . 


(64) 


Taking  the  preferred  reference  frame  Z,  to  be  a  perfectly 
known  frame  and  the  host  vehicle  body  frame  to  be  uncertain  (Za 


has  no  nominal  counterpart  but  Lt  does,  designated  here  as  X.,) 
we  have 

rjt)  =£dt)h.  (65) 


Even  though  2^  is  uncertain,  _rb  ,  the  location  of  the  antenna  phase 
center  in  the  host  vehicle  frame,  is  taken  to  be  known  perfectly. 

We  then  have 

8r4(U-  r„(t)-  £-(t)rJ  -  CMt)r£  .  (66) 


We  also  note  that  _r^  corresponding  to  the  nominal  trajectory  is 
also  perfectly  known  and,  in  fact,  is  equal  to  r1'  .  The  antenna 
phase  center  for  the  nominal  host  frame  has  the  same  relative 
location  as  the  actual  antenna  phase  center  for  the  actual  host 
frame.  Then  in  the  process  of  defining  the  skew-symmetric  matrix 
of  attitude  error  (the  difference  between  2^,  and  Xbri )  , 

8r  =  Ca(t)  -Ca  ( t )  1  rbn  =  Ca  (t)[  Cbn(t)  -  I  |rbn 

A  — b  An  — tn  — -r:  — ■  J- Ar. 

54*"  ( t )  a  r  ( t )  a  Cbn  ( t )  -I 
8r  (t)  =  Ca  (t)  $Pn(t)Cfcn(t)Ca  ( t )  rbr‘  ( t ) 

“A  =bn  ==  =a  =br,  — Ar' 

=  &¥a(t)ra  (t)  =  8¥(t)r  (t)  . 

in  i»i  -—An  msc  — Ar; 


Using  the  skew-symmetric  form  R  (t)  for  vector  r  (t)  ,  we  have 

as/^  —An 

from  (11) 


5r  ( t )  =  -  R  ( t )  ( t )  , 

A  =An 

o  -rA,n<t)  rAyn(t) 

R  =  rA:n(t)  0  ~rAxJt) 

=Al, 

"rAyr.(t)  rAxn(t)  0 


(68) 
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We  then  have  for  the  discrete  time  measurement  residual  process 


5z  (t;)  =  llrdtj  ||  -||r  .  (t.)  -8_r  (t,)  +  R  ( t , )  50  ( t , )  |j  +  8B(t.) 

J  -i  ■-1  “AT: 

{llr.  ( t ,  -T)  II  -  ||  _r  ( t ;  -T )  -5ru(t,-T)  +  ( t .  -T)  5&(  t  _  -T)  ||  *&B.At  -T)  } 


+  VPRj  (tt) 


69) 


We  define  3  x  n  matrices  Kr  and  K  and  row  vector  (1  x  n  matrix) 
such  that  we  can  relate  to  "?he  full  error  state  vector  at 
measurement  time  t , . 


6  r_  ( t , )  =  K  5  x  ( t , ) 

M(tJ  =  K!5x(t;)  (70) 

5Bu(t,)  «  JigSx ( t , ) 


At  measurement  time  tj-T,  we  obtain  the  corresponding  quantities 
but  using  the  backward  transition  matrix 

5ru(t;-T)  «  K  [<J>(tI-T/tI)8x(tl)  +wri(ti-T)] 

60.(t;-T)  =  KjfcU.-Tht.jSxdt,)  +  wD  ( t ;  -T)  ]  (71) 

5Bu(t;-T)  =  k‘[<I>(t1-T,t1)5x(t1)  +wp(ti-T)]f 


where  Wc ( 1 4 -T )  is  the  driven  response  at  t.-T  due  to  the  presence  of 
the  white  noise  in  (43)  during  the  backward  interval  from  t,  to 
t,-T.  ^p(t  -T)  ,  by  virtue  of  the  white  noise  in  the  continuous 
model,  is  a  white  noise  sequence.  It  is  here  that  uncertainty  in 
the  dynamics  model  is  introduced  in  the  measurement  model. 
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The  discrete  delta-range  measurement  residual  process  is  now 


8z  ( t i )  =  II  r  (t  j )  ||  -||  r  (t.)  - 


K  -R  (t.)K  8x(t  ) 


-j|l  r  .  ( t ;  -  T)  I!  -II  ri  (ti-T)  -  K^"  Jt,-T)K^]  -T)  ,  t , )  5x  ( t  )  +w.  (t.-T) 

+kR  {  ^"^(tx-T/tjJSxtt,)  +  Wp(t,  -  T)  | 


+vdr:  (ti)  • 
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Except  for  the  noise  introduced  by  uncertainty  in  the  dynamics 
model  in  the  form  of  Wrdtj-T),  this  is  in  the  form  given  in 
Appendix  A 

5zdr,  (tj  =  hj'  [8x(t.)  ,x(t,)  ,  t.]  +  vDK;  (73) 


This  form  is  obtained  by  boosting  vDR(t,)  by  an  appropriate  amount 
and  removing  w^t.-T)  .  In  reality,  this  is  accomplished  by  an 
increase  of  the  corresponding  term  in  the  measurement  noise 
covariance  matrix.8 

From  Appendix  A  we  have 


8z 

- DP.J 


(t.)  = 


3h, 


( 5x, x, t . : 
d  (dx) 


5:<=  8x(tJ  + 


V. 


DR] 


h .  o  .  t 


!74) 


where 


9hd  ( 8x,  x,  t , 


d  ( 6x) 

observation  matrix. 


=  hj"  [  tj ,  x  ( t , )  ]  ,  is  the  j  th 


row 


of  the 


8  It  is  to  be  noted  that  w:(ti-T)  has  the  argument  t.-T.  The 
measurement  noise  covariance  matrix,  which  needs  to  be  modified, 
however,  has  the  argument  t,.  This  added  noise  contribution  does 
not  correspond  to  any  physical  noise  process  in  the  GPS  receiver 
but  reflects  additional  uncertainty  in  the  adequacy  of  the 
resulting  measurement  model  introduced  by  using  the  assumed 
dynamics  model  in  the  development  of  the  measurement  model . 
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Noting  that  II  £. ,  ( 1 1 )  ||  and  ||r  (t1-T)||  are  not  functions  of  5x  and 

taking  wn(t,-T)=P  (by  accounting  for  this  additional  uncertainty  in 
measurement  noise  covariance  matrix) , 


IlT "  [ti;x  (tj 


3 II  r  { t  i )  -  [K  -R  (t,)K  ]  Sx| 

— 3  1  *P  ■*An  •  ***  — 

3  ( bx) 


3 II  r  ( t  -T)  -  K  -  R  (ti-T)K  ]<I>(tl-T,ti)5x|| 

- JWl - 

3kni  -^(ti-T.t^lSx 

- 3TSxl - 1 • 


8x«0 


(75) 


By  using  the  results  of  Appendix  B,  and  since  for 
5x=0/x=x„  and  r\=r_ 


h™ 


rT  (t i )  [K  -  R  ( t  )  K  ] 
tt^-xjt.)  ]  =  _1_  "p  “An 


r.  (t, 

— jn 


rT  (t  -T)  [K  -  R  (t, -T)K  ]0(t. -T,  t.  ) 

— :n  1  ~p.  *An  1  «•* 


r  (t  -T)  || 

— j  n  1 


+k •  [I  -  ^(t.-T,  t;)  ] 


(76) 


As , 


r\  (t) 

HfltTI 


=  e;,;(t) 


hT"  [t,;x  (tt)  ]  =  eT  (t,)  [K  -R  (tJK  ] 

—J  1  — n  1  — jn  i  =P  «=An  ‘  =♦ 

-  ejn(t1  -T)  [K  -  R  (t,  -T)  K  JO( t j  -T,  tt) 
+  iSs[I.  ~  tj)  ]  . 


(77) 


If  we  approximate  the  backward  transition  matrix  by  the  first 
order  form 


^(t-ThtJ  =  ^  -  F[t:;x,  (tt)  ]  T  , 


(78) 
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then 


hT"  [  t  j ;  x  (t  J  ]  =  ej  ( t ; )  [K  -  R  (t,)K  ] 

— J  1  — n  1  — jn  1  »=p  *=An  1 

-eT  ( t  i  -  T)  [K  -  R  (t,-  T)K  ]  [I  -  F[ti;x  <t2)  ]  ]T  (79) 

— jn  >s=p  a^An  1  ==♦  =  =  — n 

+k^F[t1;Xn(t1)]T. 

The  3xn  matrices  K  and  K  and  the  row  vector  kT  that  extract 

■=p  ==♦  — B 

the  position  error,  attitude  error  and  user  clock  error  from  the 
full  error  state  vector  necessarily  depend  on  how  the  full  error 
state  vector  is  constructed  but  are  readily  derived. 

Rearranging  (79),  we  have  for  the  row  corresponding  to  the 
delta-range  measurement  residual  for  the  jth  satellite 

hT"  [ t  ; x  (t2)  ]  =  [eT  (t  )  -eT  (t2-  T)  ]  K 

— -j  1  — n  1  — jn  1  — “jn  1  =sp 

-[eT  (ti)R  (t,)  -eT  (tr  T)R  (trT)]K 

— an  1  1  — an  1  =An  1  ==♦ 

+eT  (tr  T)K  F  [  t : ;  x  (t,)  ]T  (80) 

—3  n  1  1  — n  1 

-e:Tn(tl-  T )  R^  ( t  j  -  T)K  F[ti;xn(ti)  ]T 
+iiBTF[t1;xn(t1)]T  . 


The  best  way  to  compare  this  method  to  the  first  method  is  to 
see  how  the  error  state  vector  is  coupled  to  a  delta-range 
measurement  residual,  by  using  £a  =  ECEF  as  an  example,  we  briefly 
examine  the  effect  of  each  of  the  five  terms  of  (80)  for  the 
product  h" J  [ t; ;  xn  ( tj )  ]  5x  ( t; )  .  We  then  compare  with  that  which  would 
be  obtain  using  a  row  of  (25)  that  corresponds  to  a  delta  range 
measurement  residual.  To  obtain  some  of  the  terms  in  (80),  of 
course,  we  must  have  a  dynamics  matrix. 

For  the  first  term  of  (80)  we  have 

[eT  (t  )-eT  (t2-  T)  ]  K  5x(t2)  =  [eT  (tl)-eT  (t2-  T)  ]  5r  (t2) 

=  AeT  ( t  i )  5r  (tj  . 

— ]  n  1  — u  1 


The  result  is  identical  to  that  for  the  corresponding  term  in  (25) . 
In  this  case  there  is  no  problem  as  to  which  time  to  evaluate  any 
part . 
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For  the  second  term  of  (80)  the  result  is  the  same  as  for  the 
corresponding  term  in  (25)  except  here  the  time  of  evaluation  is 
quite  clear. 


The  third  term  of  (80)  requires  the  dynamics  matrix.  The 
needed  terms  in  F[tif*x  (tj)]  are  simply  understood  in  that  the  time 
derivative  of  position  errors  is  equal  to  the  velocity  errors  for 
the  frame  under  consideration.  Again,  we  have  the  same  results  as 
the  corresponding  term  in  (25)  except  here  the  time  of  evaluation 
is  given. 

e^n  ( t  i  -  T)K  F  [  t  j  ;3£n  (tj )  ]T8x(ti)  =  eJJtj-  T)  T5vu(ti)  .  (83) 


The  fourth  term  of  (80)  results  in  coupling  to  gyro  related 
errors  that  may  be  included  in  Sx.  If  all  of  these  errors  were 
lumped  into  a  single  angular  velocity  error  of  the  body  frame 
relative  to  an  inertial  frame,  it  should  be  evident  that  the 
results  for  (80)  and  (25)  would  be  the  same  except  again  for  the 
time  of  evaluation  being  explicitly  known  here.  In  this  case, 

eT  (tj-  T)R  (t,-  T)KF[t,;x  (t,)]T-)eT  (tr  T)R  (t1-T)T5tOv.  (84) 

— jn  1  =An  1  — n  1  — A  =An  1  — ^ 


The  last  term  of  (80)  yields  results  identical  to  those  in 
(25)  ,  the  part  of  the  dynamics  matrix  required  here  being  quite 
simple  in  that  the  rate  of  change  of  user  clock  offset  is  equal  to 
the  user  clock  drift. 

=T8fu.  (85) 


We  have  then  compared  the  results  of  the  second  method  to  that 
of  the  first  (for  a  delta  range  measurement  residual)  when  a  first 
order  model  is  considered.  Each  case  resulted  in  five  terms  (when 
we  lump  all  gyro  errors  into  a  single  angular  velocity  error)  and 
the  five  terms  have  been  shown  to  be  equal  except  that  the  time  of 
evaluation  is  quite  clear  for  the  second  method.  This  overcomes 
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the  major  weakness  in  the  first  anproach.  The  second  method,  as 
previously  noted,  increases  uncertainty  in  the  measurement  model 
for  any  inadequacies  in  the  dynamics  model  used.  E~,h  methods  may 
be  extended  beyond  a  first  order  model. 

The  mathematical  development  for  the  second  method  has  been 
somewhat  more  concise  than  that  for  the  initial  method.  We  have 
been  able  to  keep  track  of  time  arguments,  making  a  clear 
distinction  between  the  continuous  argument  "t"  and  discrete 
measurement  times  "t,".  In  addition,  those  quantities  evaluated 
"on  the  nominal  trajectory"  have  been  identified  and  tracked 
throughout  this  development . 
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4.0  SIMPLIFIED  ANALYSIS  OF  THE  DELTA-RANGE  MEASUREMENT  FOR  A 

S PINNING  BODY 

4.1  Introduction 

This  section  provides  a  mathematical  analysis  of  attitude 
estimation  addressing  specifically  a  spinning  body  scenario.  As  a 
result  it  facilitates  gaining  certain  insights  to  and  identifying 
sensitivities  in  spinning  body  applications.  The  analysis  is 
limited  to  the  case  of  an  initial  (constant)  attitude  error  in  a 
rotating  but  not  translating  vehicle.  The  results  reveal  the 
sensitivity  to  spin  rate,  DR  integration  time,  lever  arm  length, 
and  orientation  of  the  satellite.  The  results  also  include  an 
upper  bound  of  the  accuracy  with  which  attitude  errors  can  be 
estimated . 


4.2  Analysis  Of  Delta-Range  Measurements 


Neglecting  the  effects  of  clock  bias,  atmospheric  delays,  and 
noise,  the  doppler  shift  of  the  incoming  GPS  signal  is, 
approximately. 


w. 


DOP. 


(t; 


2Jtf  RF  d 
c  "eft 


(86) 


where  fRF  is  the  GPS  signal  carrier  frequency,  jt,  is  the  Line  of 
Sight  (LOS)  vector  (see  figure  1),  c  is  the  speed  of  light,  and  the 
double  lines  denote  magnitude  of  a  vector  quantity. 

The  delta-range  measurement  is  physically  obtained  by  sensing 
and  integrating  the  doppler  shift  in  (86)  over  an  interval  of  time 
T. 


JDR. 


4<ti>  = 


27tf 


RF 


/ 

t.-T 


^DOP-  ^ 


(87) 


Using  (86)  and  the  fundamental  theorem  of  calculus  in  (87)  we 
obtain 


zdr  (t,)  = 


/ 


-ll!  r  ( x )  |  dx 
c  dx  ~ J 


■RF  t  _T 

=  |  r .  ( 1 1 )  ||  -  ||  r^  (t,-T)  ||  . 


(88! 


Without  loss  of  generality  we  can  assume  an  inertial  frame  XYZ 
as  shown  in  figure  2.  For  simplicity  we  assume  no  motion  other 
than  a  constant  spin  (0  about  the  Y-axis.  Suppose  the  antenna  is 
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Satellite 


initially  located  on  the  X-axis.  Then,  the  lever-arm  vector  rt  is 
represented  in  the  XYZ  frame  by 

L  COS  ((Ot  :  ) 

_r(t,)=  0  <  (89) 

L  sin  (0)t  j ) 

where  L  is  the  magnitude  (length)  of  r A.  The  Line  of  Sight  (LOS) 
range  vector  from  the  j-th  satellite  to  the  antenna  at  time  t.  is 
given  by  (see  Figure  1) 

Ijft.)  =  r^.ttj)  +  rA(tj).  (90) 


In  practice,  T  is  sufficiently  small  so  that  in  the  interval 
(tj,  t j-T)  the  vector  can  be  considered  a  constant,  i.e.,  we  can 
neglect  the  effect  of  satellite  motion. 

For  the  purpose  of  this  analysis,  let  us  consider  the  LOS 
vector  as  a  function  of  the  continuous  time  variable  t  as  opposed 
to  the  discrete  time  tj.  Then 

I^UHI  =  [(rcjx+Lcos<at)2  +  ( rCjy~0 ) 2  (91) 

+  (rC3,+Ls incot ) :] 1/2  . 


Then 


[( rcjx+Lcoscot )( -Lcos incot )  +  ( rCJ.+Ls incot )  (Lcocoscot )  ] 


rj 


(92) 


Typically,  the  lever  arm-length  L  is  small  compared  to 
distance  between  the  vehicle  and  a  GPS  satellite.  Then 


rcjx  |  >>  |  L  coscot 


(93: 


rc]zl  >>  I  L  s incot  |  . 


:94) 


Since  L  is  assumed  to  be  small,  by  taking 

«r:(t)«  -  ||rc]|| 


(95) 
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for  any  time  t  and  neglecting  terms  in  Lr  ,  (92)  reduces  to 

d  „  „  -r ■  Lcosincot  +  r..,,Lcocos(Ot 

^||r,,t)|l  -  - T_r_^ - 

=  Leo  [  -e_xsincot  +  e,, ..coscot]  , 


(96) 


where  ecjx  and  ecj,  are  the  X  and  Z  components  of  a  unit  vector  e,:  in 
the  direction  of  rtj .  Then  from  (86) 

wi»rj(t)  =  —  —  Leo  [-ecjxsincot  +  eCJ_cos(ot].  (97) 


From  (97),  we  conclude  that  the  doppler  shift  is  sinusoidal  in 
time.  The  doppler  shift  zero  crossings  occur  whenever 


tan  (tot ) 


(98) 


i.e.,  whenever  the  antenna  crosses  the  projection  of  r\;j  onto  the 
X-Z  plane  (indicated  as  line  aa '  in  figure  2) . 


From  (87)  and  (97)  we  obtain 

COT 


COT 


:0Rj(t)  =  2L  sin(_^_)  [-ecjxsin(tot-_2_ 


COT 


+  e^cos  (tot-^i)  ]  .  (99) 


(99)  can  also  be  written  in  terms  of  the  inner  product  between 
the  vectors  and  u  where 


u 


-sin  (cot-i^E ) 


0 


cos  (cot  - 


COT  , 

~T ' 


(100) 


The  vector  u_  is  orthogonal  to  the  vector  denoting  the  antenna 
location  at  the  midpoint  of  the  delta-range  integration;  i.e.,  u  is 
orthogonal  to 

Therefore,  for  a  constant  rotation  rate  CO,  u.  is  parallel  to 
the  vector  Ar^,  which  denotes  the  difference  in  antenna  location  at 
the  start  and  end  of  the  delta-range  integration. 
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(iYT 

Lcos  ) 


0 

Lsin  (cot  - 


COT 

~T ' 


101) 


Ar  a  ( t )  =  r  a  ( t )  -  r  A  ( t  -T)  . 


Thus  (99)  can  be  rewritten  as 


filT 

zDRj(t)  =  2L  sin(^i)  cos  (a[ecj/ArA(t)  ] )  , 


102) 


(103) 


where  oc[e^j(  ArA(t)]  is  the  angle  between  the  unit  vector  e,3  and 
Ar,A-  As  indicated,  this  angle  is  a  function  of  e and  Ar^  ( t )  . 

The  following  observations  can  now  be  made: 

1)  From  (99),  the  delta-range  measurement  is  sinusoidal  in 
time.  The  delta-range  measurement  zero  crossings  occur  whenever 

tan(co(t--I)  )  =— 121  ,  (104) 

^  ecjx 


i.e.,  whenever  the  angle  swept  by  the  antenna  lever  arm  during  the 
integration  time  T  is  symmetrically  distributed  about  the  line  aa'  . 
This  condition  is  illustrated  in  figure  3a. 

2)  Setting  the  time  derivative  of  (99)  equal  to  zero  we 
conclude  that  the  peak  values  of  the  delta-range  measurement  occur 
whenever 


tan  (a>(t-.£ )  ) 


(105) 


i.e.,  whenever  the  angle  swept  by  the  antenna  lever  arm  during  the 
integration  time  T  is  symmetrically  distributed  about  the  line  bb' 
which  is  orthogonal  to  the  projection  of  r\;  on  the  X-Z  plane.  This 
condition  is  illustrated  in  figure  3b. 
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Figure  3 


a)  Zero  DR  measurement  and  peak  attitude  error 
observability 


b)  Peak  DR  measurement  and  zero  attitude  error 
observability 


.  Illustration  of  conditions  for  zero  and  peak  delta-range 
measurement  and  attitude  error  observability. 
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3)  From  (103)  we  conclude  that  the  maximum  peak  value  of  the 
delta-range  measurement  is  2L.  This  value  is  obtained  when  the 
vectors  e,,  and  An,  are  parallel  and  the  antenna  rotates  180  degrees 
starting  and  ending  on  the  line  aa ' .  On  the  other  hand,  the  delta- 
range  measurement  is  zero  if  the  vectors  e_.,  and  Ar*  are  orthogonal. 

4.3  Analysis  Of  Attitude  Error  Observability 

In  this  paragraph  we  investigate  the  observability  afforded  to 
attitude  errors  by  the  delta-range  measurement.  Attitude  errors 
are  directly  observable  whenever  the  corresponding  entries  in  the 
observation  matrix  H  are  non  zero  -  i.e.,  the  entries  in  column  8<J> 
of  (32).  For  "good"  observability,  these  entries  should  be  as 
large  as  possible  so  that  the  contribution  of  attitude  errors  to 
the  DR  residual  is  significant. 


We  now  return  to  a  discrete  time  formulation.  From  (32),  the 
attitude  error  8j£  is  related  to  the  delta-range  residual  for  the  j- 
th  satellite  via 


-  (AeT  R  +  eT  AR  ) . 

-“Jn«An  3 n  =A.n 


(106) 


For  our  choice  of  reference  frame  and  short  periods  of  time 
the  change  in  the  unit  vector  e.3n  is  small  so,  we  can  assume  that 
Ae^  is  approximately  zero.  Then  (106)  reduces  to 


-eT  AR 
— m  • 


(107) 


The  contribution  of  the  attitude  error  _5^_  to  the  DR 
measurement  residual  is  obtained  by  multiplying  both  sides  of  (106) 
by  8^.. 

H  ,,84)  =  -eT  AR  86  (108) 

— jOf— <-  — jn  mmpn *- 


-eT  AR  8b  -e7n 

— J>‘  =An — 3 


—  T  .  /  r  — 


Ur.(t, 


-  r. 


5<p ) 


(io9: 


or 


ej  AR_  8£<->  -e5n-  i  [ArA(T)  ]  x  80  )  . 


(110) 


From  (110),  we  conclude  that  the  contribution  of  the  attitude 
error  to  the  delta-range  residual  is  maximum  when  the  vectors 
e.:n,  Ar^,  and  are  mutually  orthogonal.  The  contribution  is  zero 
if  any  two  are  in  the  same  plane. 
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From  figure  2  we  conclude  that  the  skew-symmetric  matrix  RA. 


is 


0  -Ls  incot,  0 
R  a  Ls  incot,  0  -LcosCOt, 
0  LcosCOt,  0 


(111) 


Then 


AR  = 


0  -L  (sincot1-sinco(t1-T)  )  0 

-sym  0  -L  ( coscot  1  -cosco(  t .  -T)  ) 

0  -sym  0 


112) 


where  "sym"  denotes  an  entry  equal  to  the  symmetric  entry  in  the 
matrix.  Using  trigonometric  identities  (112)  reduces  to 


AR  = 

“An 


fiVT  T 

0  -2Lsm_^icosco(  t , -.j  ) 


-sym 

0 


0 

-sym 


cor  T 

2Lsin___sinco(  t,  -_ , 

^7  ■  T 


0 


.  (113) 


Since  we  have  assumed  that  L  is  small  compared  to  ,  e.:  is 
approximately  equal  to  e j .  Finally,  assuming  we  are  only  interested 
in  the  y-component  of  the  attitude  error,  we  obtain  from  (107), 


fiYT  T 

h:6^  =  2L  sin(_^i)  [ecjxcosco(t,-^.] 


+  ec,.sinco(t,--y)  ]  .  (114! 


From  (114)  we  conclude  that  h_Sf(.  is  zero  if  any  of  the  following 
three  conditions  occur: 


L  =  0  , 


(115) 


COT 

~r 


NJC 


(N  =  0,1,2,  .  .  .), 


(116) 


ecjxcos  (cot 


COT 

'~~T 


+  e_,.sin  (cot , 


COT 

~T 


=  0 


(117) 
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Obviously,  the  condition  in  (115)  can  occur  if  and  only  if  the 
lever  arm  is  zero.  The  condition  in  (116)  can  occur  if  there  is  no 
rotation  or  if  the  rotation  rate  is  such  that  a  multiple  of  360 
degrees  is  swept  during  the  delta-range  integration  time  T. 

The  condition  in  (117)  can  occur  if 

e  .  =  e...  =  0  (118) 

..  .  f. 

or 

tan  (tot ,  -i£E  )  =  -—1 1  .  (119) 

^  ec3  = 


The  condition  in  (118)  occurs  if  the  satellite  is  directly 
above  (or  below)  the  vehicle;  i.e.,  if  the  vector  ru,  in  figure  2 
has  only  a  Y  component .  The  condition  in  (119)  occurs  if  the  angle 
swept  by  the  antenna  lever  arm  during  the  delta-range  integration 
interval  is  symmetrically  distributed  about  the  line  bb'  in  the  X-Z 
plane  -  see  Figure  3 . 

We  can  express  (114)  in  terms  of  the  inner  product  between  eCJ 
and  the  antenna  location  vector  at  the  midpoint  of  the  delta-range 
integration,  rA(ti-T/2). 


h 


=  2  sin(i^E)  [eh  rA(tr^)  ] 

=  2  L  sin(i^E)cos(P[ecj,rA(tJ-^)  ]  )  , 


(120) 


where  (3  [  ,  rA(t,-T/2)]  is  the  angle  between  the  unit  vector  e,.  and 

the  antenna  location  vector  rA(t,-T/2). 

From  (12  0)  we  conclude  that  hj6^y  is  sinusoidal  and  has  a 
maximum  peak  value  of  2L.  This  value  is  obtained  when  the  vectors 
e^  and  _rA(t1-T/2)  are  parallel  and  the  antenna  rotates  180  degrees 
starting  and  ending  on  line  bb'.  hjfi4y  is  zero  if  the  vectors  e... 
and  rA(ti-T/2)  are  orthogonal.  It  follows  that  maximum  attitude 
error  observability  is  obtained  when  the  DR  measurement  is  at  a 
zero  crossing. 

4.4  Attitude  Estimation  Error  Bound 


Suppose  an  Information  Kalman  Filter  is  used  to  estimate  a 
constant  attitude  error  8<{>y  using  DR  measurements.  The  variance  of 
the  estimation  error  at  tk  is  (Anderson  &  Moore,  1979) 
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(  121) 


(°  c6*  (  t  ,)  )  =  (0-£ 


h  6*  (t,)‘ 


where  oDR  denotes  the  sigma  value  of  the  observation  noise. 

Numerical  values  for  h,8^y(t1)  are  required  in  order  to  evaluate 
the  sum  in  (121)  .  Assuming  that  tk  is  sufficiently  large,  we  can 

assume  that  the  sum  can  be  evaluated  using  an  average  (and 
therefore  constant)  value  in  place  of  h }  8+y  ( 1 4 )  ;  i.e., 


£  Vit,) 


-  (k  +  1) 


fAvg(hjS^)); 


(122) 


Assuming  that  to  a  sufficient  approximation 


Avg(hj84 


Peak  (hj8>  ) 


:123) 


(121)  reduces  to 


«T2e8*  (tk))':  =  (t0))-5  -  (k+1) 


Peak  (h,64  )  - 


*  ( 


,  Peak  (h-5#  ) 2 

*  — ..  ■ 


>  Peak<hi8t,); 


(124: 


Using  (120)  in  (124)  and  taking  the  square  root  of  both  sides 
we  obtain 


<Wtk)  - 


L  sin  ( i!^  )  cos(3 


;125) 


for  t„>>0 
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The  following  conclusions  can  be  drawn: 

1)  The  uncertainty  in  estimating  attitude  errors  0e8*  with  an 
offset  antenna  is  proportional  to  the  delta-range  measurement 
noise,  i.e.,  oDR.  (0DR  is  inversely  proportional  to  the  SNR  at  the 
input  of  the  carrier  tracking  loop) . 

2)  0e8<(  is  inversely  proportional  to  the  length  of  the 
antenna  lever  arm  L. 

3)  0e8i  is  minimized  if  the  angle  swept  during  the  delta- 
range  integration  is  180  degrees. 

4)  0e8*  is  minimized  if  the  satellite  is  in  the  plane  defined 
by  the  rotating  antenna  and  the  angle  swept  during  the  delta-range 
integration  is  symmetrically  distributed  about  the  vector  rcj . 

We  note  that  (125)  is  only  a  bound.  As  a  result,  it  is  not 
very  useful  in  determining  the  uncertainty  in  the  attitude 
estimates  exactly.  It  is,  however,  useful  in  evaluating  the 
sensitivity  of  the  uncertainty  to  L  and  COT.  To  illustrate  this, 
let  7  be  a  factor  such  that  strict  equality  is  achieved  in  (125) 
for  some  L0,  co0,  and  T0;  i.e., 


k ;  L0 ,  G)0 ,  T0 )  - 


7  o 


DR 


C0.T 


L0  sin(  )  cosp. 


(126) 


Then,  for  different  L,  0),  and  T  we  have 


_  Lo 


sin  ( 


<*>oT0 


0e8,  ( t  k ;  L ,  CO,  T)  = 


L  •  ,  COT  . 

sin  ( ) 


°c6*  ( t  k ;  L0 ,  C00 ,  T0 ) 


(127) 


For  example,  increasing  the  lever  arm  from  L  =  5  inches  to  20 
inches,  reduces  the  steady  state  attitude  error  bound  by  a  factor 
of  4.  On  the  other  hand,  with  0%  =  to  =  450  deg/sec,  increasing  the 
DR  integration  time  from  T  =  0.78  sec  to  T  =  1.0  sec,  reduces  the 
steady  state  attitude  error  bound  by  a  factor  of  9. 
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5.0 


SIMULATION  RESULTS 


5.1  Introduction 

This  section  provides  covariance  simulation  analysis  for  a 
spinning  body  by  using  the  more  general  results  of  section  3.  The 
simulation  scenario  consists  of  a  vehicle  falling  freely  in  an 
exoatmospheric  trajectory.  The  results  illustrate  the  ability  of 
estimating  attitude  and  attitude  rate  errors  with  an  offset 
antenna.  The  sensitivity  to  the  lever  arm  length  and  the  spin  rate 
were  also  investigated.  To  an  approximation,  the  scenario  here 
resembles  the  conditions  assumed  in  the  analysis  of  section  4. 
Consequently,  the  insights  gained  in  section  4  carry  over  to  the 
more  realistic  conditions  of  this  section. 

5.2  Scenario  Description 

We  consider  a  vehicle  spinning  about  its  Yb  (body  frame)  axis 
at  a  constant  roll  rate  of  co  deg/sec.  The  vehicle  is  assumed  to  be 
equipped  with  a  tightly  integrated  GPS/Inertial  navigator  mounted 
along  the  Yb  axis  and  near  the  vehicle's  center  of  gravity.  The 
GPS  antenna  is  mounted  on  the  perimeter  of  the  vehicle  at  a 
distance  L  from  the  Yb  axis .  The  antenna  therefore  rotates  at  the 
vehicle  spin  rate  co.  The  vehicle  moves  along  an  exoatmospheric 
trajectory.  The  vehicle's  altitude  profile  versus  time  is  shown  in 
figure  4. 

Assuming  a  Yb  gyro  scale  factor  error  of  50  ppm,  a  spin  rate 
of  450  deg/sec,  and  a  scenario  duration  of  1,400  sec,  an  attitude 
error  about  the  Yb  axis  of  31.5  degrees  is  expected.  The  main 
objective  of  the  simulations  is  to  demonstrate  the  ability  to 
estimate  this  attitude  error  by  using  the  approach  discussed  in 
this  report . 

5.3  Simulation  Program  Description 

The  simulations  were  performed  with  Galaxy's  Covariance 
Analysis  Package  (CAP)  which  is  shown  in  figure  5.  CAP  includes  a 
trajectory  generator  program  (TRAGEN) ,  a  covariance  analysis 
program  (GGPSIM) ,  a  sensitivity  analysis  program  (SENSIS),  and 
various  data  reduction  and  plotting  software  packages.  A  short 
description  of  these  programs  is  included  in  the  following 
paragraphs . 
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Figure  5.  Covariance  analysis  package  (CAP)  block  diagram 


5.3.1  Trajectory  Generator  (TRAGEN)  Description 

TRAGEN  is  a  6  Degree  of  Freedom  (DOF)  ,  free  falling  body 
trajectory  generator. 

5.3.2  Covariance  Simulator  (GGPSIM)  Description 

GGPSIM  is  a  covariance  simulation  program  of  a  tightly 
integrated  GPS/Inertial  navigatoi .  A  top  level  block  diagram  is 
shown  in  figure  6. 

GGPSIM  contains  a  "truth"  and  a  "filter"  covariance  model. 
The  truth  model  represents  all  errors  that  would  affect  the 
navigator.  The  truth  model  consists  of  the  80  errors  listed  in 
table  2.  Collectively,  these  errors  comprise  the  true  error  state 
of  the  navigator.  The  "filter"  model  consists  of  a  subset  of  the 
errors  in  the  true  error  state.  These  errors  comprise  the 
navigator's  Kalman  Filter  state  vector,  i.e.,  they  represent  the 
system's  perception  of  the  errors  in  the  real  world.  The  errors 
that  are  typically  included  in  the  filter  state  vector  are  the  ones 
that  can  be  unambiguously  estimated  (i.e.,  ai a  fully  observable) 
through  GPS  measurements.  The  maximum  number  of  such  states  in 
GGPSIM  is  56  and  are  listed  in  table  2.  In  practice,  the  number  of 
errors  included  in  the  filter  state  vector  is  further  limited  due 
to  processing  constraints.  GGPSIM  allows  the  user  to  arbitrarily 
select  the  states  in  the  filter.  This  feature  allows  analysis  of 
any  conceivable  GGP  Kalman  Filter  configuration 

GGPSIM  is  initialized  with  the  covariance  of  the  80-element 
true  error  state  and  the  covariance  of  the  filter  state.  These 
covariances  constitute  the  true  ana  GGP  Kalman  Filter  "error 
budget".  The  truth  and  filter  covariances  are  propagated  forward 
in  time  according  to  the  vehicle's  trajectory  until  a  set  of  GPS 
measurements  is  obtained.  GGPSIM  simulates  a  maximum  of  10  GPS 
channels,  supplying  the  filter  with  9  pairs  of  pseudorange  and 
delta-range  measurements.  (It  is  assumed  that  one  channel  is  used 
for  simultaneous  ionospheric  delay  measurements,  and  to  aid  in 
acquisition.)  Therefore,  for  the  10-channel  receiver,  only  nine 
channels  provide  pseudorange  and  delta  range  measurements) .  The  10 
satellites  are  selected  according  to  first  come  first  served  from 
all  in  view.  GGPSIM  also  allows  simulation  with  a  5  channel  GPS 
receiver  (4  pairs  of  pseudorange  and  delta  range  measurements) .  In 
this  case  the  satellites  are  selected  according  to  GDOP.  The 
measurements  are  computed  from  a  24  GPS  Satellite  Constellation 
model  which  is  integrated  with  GGPSIM. 

The  measurements  are  processed  with  a  Kalman  Filter  which 
updates  the  filter  covariance  at  1  Hz  with  Carlson's  algorithm 
(Anderson  &  Moore,  1979;  Pitman,  1962) .  The  computed  Kalman  Filter 
gains  a^e  used  in  configuring  (through  augmentation)  the  gain  for 
the  truth  model. 
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Figure  6.  GGPSIM  block  diagram. 
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Table  2.  GGPSIM  true  and  filter  error  states 


STATE  DESCRIPTION 


Position  Errors  _ 


Velocity  Errors 


IMU  Attitude  Errors 


Gravity  Deflections  and  Anomalies _ 


Gyro  Drift  Rates 


Gyro  Input  Axis  G-Sensitivities 


Gyro  Spin  Axis  G-Sensitivities 


Gyro  Output-Axis  Ga  Sensitivities 


Gyro  Scale  Factor  Errors 


Gyro  Input  Axes  Misalignments 


Accelerometer  Biases 


Accelerometer  Scale  Factor  Errors 


Accelerometer  Input  Axes  Misalignments 


Altitude  Sensor  Error 


Clock  Errors: 

Phase  Error 
Frequency  Error 
Aging  Error 

Random  Frequency  Error 
Acceleration  Sensitivities 


Satellite  Ephemeris  and  Clock  Residual  Ranging 
Errors 


TOTAL  NUMBER  OF  STATES 


TRUTH 

MODEL 


FILTERS) 


(*)  Maximum  Number  of  States 
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The  truth  model's  covariance  is  updated  with  the  Joseph  form.  The 
covariances  are  then  propagated  forward  in  time  using  trajectory- 
data  until  the  next  set  of  GPS  measurements  is  obtained. 

Two  GGPSIM  features  are  worth  noting.  First,  the  filter  gains 
are  computed  using  only  the  filter  covariance  elements.  This 
feature  enables  GGPSIM  to  simulate  realizable  Kalman  Filters. 
Second,  GGPSIM  automatically  adjusts  the  filter  input  noise  process 
covariance  (typically  denoted  by  Q)  to  account  for  the  true  error 
states  that  were  not  included  in  the  filter  state.  This  feature 
enables  GGPSIM  to  simulate  fine-tuned  filters. 

5.3.3  Sensitivity  Analysis  Program  (SENSIS) 

SENSIS  is  a  simulation  data  post  processing  program  that 
identifies  the  sensitivity  and  correlation  of  navigation  errors 
(e.g.,  position  errors,  velocity  errors,  etc.)  to  inertial  system 
errors  (e.g.,  initial  errors,  alignment  errors,  biases,  noise 
processes,  etc.).  In  this  analysis,  SENSIS  was  used  to  identify 
correlations  of  attitude  and  attitude  rate  errors  to  other 
navigation  errors. 

5.4  Simulation  Conditions  And  Assumptions 

This  paragraph  describes  the  conditions,  parameter  values,  and 
assumptions  used  in  the  simulations.  The  error  budget  for  the  GPS 
and  the  inertial  errors  is  shown  in  table  3.  These  values  are 
typical  of  advanced  GPS  receiver  designs  and  navigation  grade 
inertial  components. 

A  short  description  and  rationale  for  these  conditions  is 
given  in  the  following  paragraphs. 

5.4.1  I MU  Error  Budget 

A  strapdown  IMU  with  Fiber  Optic  Gyros  (FOG)  was  assumed. 
Since  FOGs  exhibit  very  little  g  or  g2  sensitivity,  these  errors 
were  assumed  to  be  zero.  Gyro  and  accelerometer  biases  and  scale 
factors  were  chosen  according  to  the  requirements  for  the  Phase  1 
GPS  Guidance  Package  (GGP) .  These  values  imply  an  inertial  system 
of  roughly  0.58  nmi/hr.  The  remaining  IMU  errors  were  assigned 
values  that  are  typical  of  this  class  of  inertial  navigation 
system. 

It  is  noted  that  the  values  assumed  for  the  inertial 
components  are  not  typical  of  inertial  systems  intended  for  space 
applications;  for  example,  space  applications  generally  require 
gyro  scale  factor  errors  of  a  few  ppm  (Pitman,  1962)  as  opposed  to 
50  ppm  assumed  here.  The  50  ppm  value  was  intentionally  assumed 
here  in  order  to  illustrate  the  concept. 
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Table  3.  Simulation  error  budget. 
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Table  3.  Simulation  error  budget.  (Continued) 
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X  92  ppm  N/A 

Accelerometer  Scale  Factors  Y  51  ppm  N/A 
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5.4.2  GPS  Error  Budget 

A  10-channel  GPS  receiver  was  assumed  to  establish  an  upper 
bound  of  performance  capability.  It  was  assumed  that  9  channels 
tracked  satellites  (all  in  view)  supplying  the  fil  r  with  9  pairs 
of  pseudorange  and  delta  range  measurements.  The  'O'?'  channel  was 
assumed  to  be  providing  ionospheric  measurements. 

The  GPS  receiver  accuracies  were  chosen  according  to  SS-US-200 
(1979)  and  SS-GPS-300B  (1980)  and  are  similar  to  those  of  the  GPS 
UE  3A  receiver.  These  values  were  chosen  as  a  worst  case  in 
accordance  with  the  minimum  requirements  specified  in  the  GGP 
System  Specif ication . 

The  clock  frequency  accuracy  was  derived  from  SS-US-200 
(1979)  .  The  remaining  clock  errors  were  assigned  values  which  are 
typical  of  oscillators  used  in  current  GPS  applications. 

The  following  assumptions  were  made:  a  mask  angle  of  5°  above 
the  horizon;  3-db  noise  figure;  omni  directional  antenna;  antenna 
shadowing  was  neglected  as  a  simplifying  assumption;  satellite 
ephemeris  and  clock  residual  errors  expected  of  the  Block  II 
constellation  were  used. 

As  previously  noted,  the  rotating  antenna  introduces  a 
systematic  phase  shift  to  the  received  GPS  signal.  It  is  assumed 
that  this  phase  shift  is  properly  compensated  for  by  the  tightly 
integrated  navigator.  This  is  a  reasonable  assumption  since  the 
gyros  are  able  to  sense  the  rotation  to  within  50  ppm.  The 
residual  phase  error  is  assumed  to  be  is  negligible  as  far  as  the 
GPS  pseudo-range  and  delta-range  measurements  are  concerned. 

5.4.3  Kalman  Filter  Configuration 

An  "optimal"  Kalman  Filter  configuration  was  used.  The  state 
vector  included  all  filter  states  listed  in  table  2  except  for  gyro 
g  and  gz  sensitivities  (a  total  of  47  states) .  This  configuration 
was  used  in  order  to  identify  which  and  to  what  extent  the  various 
error  states  can  be  estimated.  In  reality,  reduced  order  filters 
will  be  used.  In  these  realistic  configurations,  states  that  can 
not  be  estimated  to  a  significant  degree  are  removed  from  the  state 
vector.  This  analysis,  however,  was  not  done  here. 

5.5  Simulation  Results 

Results  were  first  obtained  for  three  sets  of  lever  arm  length 
L,  spin  rate  to,  and  delta-range  measurement  integration  time  T. 
These  results  are  discussed  as  Cases  1,  2,  and  3  below.  Additional 
results  were  then  obtained  to  reveal  the  sensitivity  of  attitude 
errors  to  antenna-lever-arm  length  and  spin  rate. 


54 


5.5.1  Case  1:  L  =  0.0  inches,  CO  =  450  deg/sec,  T  =  0.78  sec 

This  case  was  considered  in  order  to  establish  a  baseline. 
Physically,  this  case  would  result  if  the  antenna  were  mounted  on 
the  Yb  axis  -  figure  1.  As  per  the  discussion  in  section  2, 
attitude  error  estimation  with  a  zero  lever  arm  is  not  possible. 

Figure  7  shows  the  resulting  Position  Errors  vs  Time  and 
figure  8  shows  the  resulting  Velocity  Errors  vs  Time.  The  results 
indicate  good  estimation  of  position  and  velocity  errors  due  to  the 
processing  of  GPS  measurements.  The  residual  position  errors  are 
due  to  satellite  biases.  The  residual  velocity  errors  are  small 
(0.01  m/s)  due  to  smooth  dynamics,  good  satellite  observability, 
and  lack  of  jamming.  The  "bumps"  in  the  velocity  error  profiles 
are  due  to  satellite  switching. 

In  typical  terrestrial  navigation  applications,  good  position 
and  velocity-error  estimation  would  imply  good  attitude  error 
estimation  since  attitude  errors  are  coupled  to  velocity  errors  via 
the  specific  forces  (Schuler  Loop) .  In  this  case,  however,  the 
vehicle  is  under  free-fall  conditions  and  the  specific  forces  are 
zero.  Consequently,  the  attitude  errors  are  decoupled  from 
position  and  velocity  errors  and  without  a  lever  arm  there  is  no 
mechanism  for  attitude  estimation. 

Due  to  the  scale  factor  error,  the  attitude  error  about  the  Yb 
axis  is  expected  to  grow.  Indeed,  the  attitude  error  profiles  in 
figure  9  (Geographic  Coordinates)  and  figure  10  (RSS  Attitude 
Error)  show  a  roughly  linear  growth.  The  attitude  error  is  shown 
to  grow  to  approximately  30  degrees  as  predicted  earlier. 

In  addition,  none  of  the  instrumentation  errors  affecting 
attitude  (e.g.,  gyro  scale  factor  errors,  gyro  bias,  etc.)  can  be 
estimated.  For  example,  figure  11  shows  all  gyro  scale  factors 
errors  remain  constant  throughout  the  trajectory. 

In  essence,  under  the  conditions  of  this  case,  GPS  succeeds  in 
estimating  the  translational  errors  (position  and  velocity)  as  if 
the  vehicle  were  a  point  mass.  There  is  no  mechanism  to  estimate 
attitude  errors  either  as  direct  observables  or  through  coupling  to 
other  directly  observable  errors. 
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Figure  7.  Case  1:  position  errors  vs.  time. 
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Figure  11,  Case  1:  gyro  scale  factor  errors  vs.  time. 
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5.5.2  Case  2:  L  =  7.0  inches,  (0  =  450  deg/sec,  T  =  0.78  sec 

In  this  case,  a  lever  arm  of  L  =  7  inches  as  shown  in  figure 
1  is  assumed.  Again,  good  position  and  velocity  error  performance 
is  obtained  (figures  12  and  13).  In  addition,  however,  attitude 
errors  are  estimated.  Figures  14  and  15  illustrate  the  attitude 
error  profiles  in  Geographic  Coordinates  and  RSS  value, 
respectively . 

The  results  show  that  attitude  error  has  been  estimated  to 
within  2  degrees  -  roughly,  an  order  of  magnitude  improvement  over 
Case  1.  The  RSS  attitude  error  profile  in  figure  15  consists  of  a 
period  of  transient  performance  followed  by  a  nearly  flat  behavior. 
This  is  due  to  the  estimation  of  the  instrumentation  errors. 
Figure  16  shows  that  the  Yb  axis  gyro  scale  factor  error  is 
estimated  to  within  TBD  ppm.  The  other  scale  factor  errors  are  not 
estimated  since  there  is  no  rotation  about  those  axes. 

Note  that  the  position  and  velocity  error  profiles  here  are 
identical  to  those  in  Case  1.  This  illustrates  that  position  and 
velocity  errors  are  fully  decoupled  from  attitude  errors  under  free 
fall  conditions.  Consequently,  although  here  attitude  errors  were 
reduced  through  estimation,  the  position  and  velocity  errors  were 
not  reduced  any  further. 
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Figure  12.  Case  2:  position  errors  vs.  time. 
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Figure  13.  Case  2:  velocity  errors  vs.  time. 
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Figure  15.  Case  2:  RSS  attitude  errors  vs.  time. 
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5.5.3  Case  3:  L  =  7.0  inches,  0)  =  450  deg/sec,  T  =  1.00  sec 

In  this  case  the  delta-range  integration  time  is  increased  to 
1  sec.  As  per  the  discussion  in  Section  4,  the  increase  in  delta- 
range  integration  time  results  in  a  more  favorable  product  coT 
(angle  swept  during  the  DR  integration  time) .  Again,  the  position 
and  velocity  error  profiles  were  idehtical  to  the  ones  obtained  in 
the  previous  cases  (figures  17  and  18)  . 

The  attitude  error  profiles  are  shown  in  figures  19  and  20  in 
Geographic  Coordinates  and  RSS  value,  respectively.  Figure  20  shows 
that  the  attitude  error  is  now  estimated  to  within  0.2  degrees, 
roughly  two  orders  of  magnitude  improvement  over  Case  1  and  an 
order  Qf  magnitude  improvement  over  Case  2.  Figure  21  shows  that 
the  scale  factor  error  is  now  estimated  to  within  0.2  ppm. 

5.5.4  Sensitivity  Results 

Additional  simulations  were  performed  using  the  lever  arm 
length  L  and  the  spin  rate  CD  as  simulation  parameters.  Figure  22 
shows  the  RSS  attitude  error  at  the  end  of  the  simulation  vs  the 
lever  arm  length.  Figure  23  shows  the  RSS  attitude  error  at  the  end 
of  the  simulation  vs  the  spin  rate.  These  results  agree  with  the 
predictions  in  Section  4,  i.e.,  the  attitude  error  is  inversely 
proportional  to  the  lever-arm  length  and  shows  a  cosecant 
dependence  on  the  spin  rate. 
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Figure  17.  Case  3:  position  errors  vs.  time. 


—  EAST  VELOCITY  ERROR 
—  NORTH  VELOCITY  ERROR 

M/S  —  VERTICAL  VELOCITY  ERROR 


Figure  18.  Case  3:  velocity  errors  vs.  time. 
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Figure  23.  Final  RSS  attitude  sensitivity  to  spin  rate 
(T  =  0.78  sec,  L  -  7  inch). 
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6 . 0  SUMMARY 


A  method  was  presented  which  permits  estimation  of  attitude 
and  attitude  rate  errors  with  an  offset  GPS  antenna  mounted  on  a 
platform  undergoing  attitude  changes.  The  method  exploits  the 
attitude  information  inherently  present  in  the  pseudo  range  and  the 
delta  range  measurements.  The  attitude  information  is  primarily 
recovered  through  processing  the  delta  range  measurement  with 
standard  recursive  estimation  techniques  such  as  the  Kalman  Filter 
algorithms  of  GPS/INS  systems.  A  brief  comparison  to  GPS 
interferometric  techniques  for  attitude  estimation  was  included. 

Certain  design  considerations  and  performance  sensitivities 
were  identified  and  verified  with  covariance  simulations. 
Specif ically ,  attitude  estimation  performance  is  sensitive  to  the 
angle  coT  angle  swept  during  the  DR  integration,  the  length  of  the 
antenna  lever  arm,  and  satellite  directions. 

The  method  is  ideally  suited  to  spinning  vehicles  in 
exoatmospheric  trajectories.  In  these  applications  recovery  of 
attitude  errors  with  conventional  GPS/INS  algorithms  is  not 
possible  due  to  the  decoupling  between  velocity  and  attitude 
errors.  Certain  benefits  may  be  possible  in  terrestrial 
applications  where  rotary  motion  of  the  antenna  platform  is 
inherent,  e.g.,  a  helicopter  with  a  GPS  antenna  mounted  on  the 
blade  of  the  rotary  wing. 

In  general,  an  integrated  GPS/INS  system  is  required  for 
"good"  performance.  In  that  case,  the  implementation  impact  is 
minimum  -  a  few  terms  in  the  observation  matrix  of  the  GPS/INS 
Kalman  filter.  A  stand  alone  GPS  receiver  may  suffice  in  certain 
applications  such  as  benign  or  well  known  nominal  angular  dynamics. 
In  that  case,  the  implementation  impact  also  includes  constructing 
an  attitude  solution  in  near  real  time. 

Although  use  of  multiple  offset  antennas  (and  their  associated 
receivers)  has  not  been  given  specific  treatment,  the  corresponding 
general  observation  matrix  is  simply  obtained  by  incorporating  a 
similar  set  of  rows  for  each  offset  antenna/receiver  combination. 
This  may  be  accomplished  whether  the  observation  matrix  is  based  on 
the  use  of  pseudo-range  or  delta-range  residuals  or  both. 
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APPENDIX  A 


Given  the  error  state  vector,  the  nonlinear  equation  for  the 
measurement  process  and  its  counterpart  for  the  nominal  trajectory 

5x  ( t )  =  x  ( t )  -  x,  ( t ) 

_z(ti)=h[x(t,),tj+v(t.)  (Al) 

Zt(t.)  =  h  (t;  )  /  t  J  • 

The  measurement  residual  is  then 

8z  { t , )  =  z  ( t , )  -zr  (t1)=h[x(t1),t.]-h  [  xn  (t,),t,]+v(t.).  (A2) 

If  we  expand  h[x(tj,  t,]  about  x^th)  using  Taylor's  Theorem 

UxUJ  ,  tj  =h[xn(ti)  ,tj  +-^l(^ti)  l^itp  U[t±)  -Xn  (tA)  ]  (A3) 

+o2  [x(tt)  -xn(ti)  ]  . 


We  then  have 

h'  [8x  (t , )  ,  x  (t . )  ,  t. ,  ]  A  h[x(;)  ,  tj  -  h[x  (t.)  ,  t. 

3h  (x,  t . ) 
cJx 


V.V  5x 


-  (5x) 


(A4) 


and,  therefore,  by  treating  8x  as  an  independent  variable 

3h'(5x,x,t1)  3h(x,t) 

~  J  (~6xi - ^  ~Jx  '  °<5*>  • 

Evaluating  at  8x  -  0 

3h'  (8x,x,  tj  ,  3h(x,t,). 

If  h[x(t,),  tj  has  the  more  restrictive  form 


(A5) 


(A6) 


h[x(tj  ,  tj  = 


h.  [x  (t, )  ,  tj 
h:.  [x(t, )  ,  tj 

|hm[x(t,)  ,  tj 


(A7) 


A- 1 


and,  similarly, 


h, '  [5x(tJ  ,x(tL)  ,  t.J 
hm '  [ 8x ( t , )  ,  x  ( t , )  ,t,] 

then  we  may  deal  only  with  a  representative  "row"  for  the  above. 
In  this  event  we  have 

Sz^tj)  =  hj  [8x(ti)  ,x(tj)  ,x,  (tj)  ]  +v.(ti) 

Vi t I*  f .v  ft.  M  -  dhj  (x,  t4)  ,  _  ahv(6x,x,  (t,)  (A9 

- 3T8xl - Ux*°  ' 


where 


Hft^-xjt,)  ]  = 


’hj”  1 1 i ; ( t , ) ] 

hj“  [ t  ; x  (t; ) ] 

hT "  [ t  ; x  ( t  j )  ] 

— m  1  — n  1 


(A10) 


_  8h,"  (5x,x,  t. )  .  ~ 

— '■  ~i  ,77 - l*x-o5i£<ti^+v’j<ti)+  h. 


o  .  t .  . 
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APPENDIX  B 


In  this  appendix,  we  show  that 


9 1|  V+Ku|| 

3jl 


VTK 

irar  ■ 


(Bl) 


We  note 
operator . 


that  — 
ou 


functions  as 


a  post 


1  ^ 

f(V  +  Ku)T  (V  +  Ku)1t _Z_ 

l  —  *=—  —  =—  J  du 

1  -  1 

-j[  ( V  +  Ku)  T  ( V  +  Ku]  T2  ( V  +  Ku)  TK 

( V  +  Ku)  TK 

"1  V""+  KuIF 

Therefore, 

ail  V  +  Ku|  VTK 
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